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Abstract 



We analyze the static QCD potential Vqcd(^) in the distance region 0.1 fm < r < 1 fm 
using perturbative QCD and operator-product expansion (OPE) as basic theoretical tools. 
We assemble theoretical developments up to date and perform a solid and accurate anal- 
^ ■ ysis. The analysis consists of 3 major steps: (I) We study large-order behavior of the 

perturbative series of Vqcd(^) analytically. Higher-order terms are estimated by large-/3o 
approximation or by renormalization group, and the renormalization scale is varied around 
the minimal-sensitivity scale. A "Coulomb" +linear potential can be identified with the 
scale-independent and renormalon-free part of the prediction and can be separated from the 
renormalon-dominating part. (II) In the frame of OPE, we define two types of renormaliza- 
tion schemes for the leading Wilson coefficient. One scheme belongs to the class of conven- 
tional factorization schemes. The other scheme belongs to a new class, which is independent 
of the factorization scale, derived from a generalization of the "Coulomb" +linear potential of 
(I). The Wilson coefficient is free from IR renormalons and IR divergences in both schemes. 
We study properties of the Wilson coefficient and of the corresponding non-perturbative con- 
tribution (5£'us( r ) hi each scheme. (Ill) We compare numerically perturbative predictions of 
the Wilson coefficient and lattice computations of Vqcd(^) when ni = 0. We confirm either 
correctness or consistency (within uncertainties) of the theoretical predictions made in (II). 
Then we perform fits to simultaneously determine 5E\js(r) and r o^-^^° P (relation between 
lattice scale and Aj^g). As for the former quantity, we improve bounds as compared to the 
previous determination; as for the latter quantity, our analysis provides a new method for its 



. 3-loop 
MS 

We elucidate the mechanism for the sensitivities and examine sources of errors in detail. 



determination. We find that (a) 5Eus(r) = is disfavored, and (b) ^0^^- = 0.574 ±0.042. 



1 Introduction 

In this article, we study the QCD potential for a static quark-antiquark (QQ) pair, in the distance 
region 0.5 GeV" 1 (0.1 fm) < r < 5 GeV" 1 (1 fm). This region is known to be relevant to the 
spectroscopy of the heavy quarkonium states. We use perturbative QCD and operator-product- 
expansion (OPE) as basic theoretical tools, taking advantage of dramatic theoretical developments 
that took place in the last decade. In addition, we use recent accurate results of lattice computa- 
tions of the QCD potential. 

For 30 years, the static QCD potential Vqcd(?") has been studied extensively for the purpose of 
elucidating the nature of the interaction between heavy quark and antiquark. Generally, Vqcd(^) 
at short-distances can be computed accurately by perturbative QCD. On the other hand, the po- 
tential shape at long-distances should be determined by non-perturbative methods, such as lattice 
simulations or phenomenological potential-model analyses; in the latter approach phenomenolog- 
ical potentials are extracted from experimental data for the heavy quarkonium spectra. 

Computations of Vqcd(?") in perturbative QCD has a long history. At tree-level, Vqcd(?") is 
merely a Coulomb potential, — (4/3)(as/r), arising from one-gluon-exchange diagram. The 1-loop 
correction (with massless internal quarks) was already computed in [1, 2]. The 1-loop correction 
due to massive internal quarks was computed in [3]. It took a rather long time before the 2- loop 
correction (with massless internal quarks) was computed in [4]; part of this result was corrected 
soon in [5]. The 2- loop correction due to massive internal quarks was computed in [6, 7, 8]; 
misprints in [7, 8] were corrected in [9]. The logarithmic correction at 3-loop originating from 
the ultrasoft scale was first pointed out in [1] and computed in [10, 11]. A renormalization-group 
(RG) improvement of VqcdM at next-to-next-to- leading logarithmic order (NNLL), including the 
ultrasoft logarithms, was performed in [12]. (There exist estimates of higher-order corrections to 
the perturbative QCD potential in various methods [13, 14, 15].)* 

For a long time, the perturbative QCD predictions of Vqcd(^) were not successful in the dis- 
tance region relevant to the bottomonium and charmonium states, 0.5 GeV" 1 < r < 5 GeV -1 . 
In fact, the perturbative series turned out to be very poorly convergent at r > 0.5 GeV -1 ; uncer- 
tainty of the series is so large that one could hardly obtain meaningful prediction in this distance 
region. Even if one tries to improve the perturbation series by certain resummation prescriptions 
(such as RG improvement), scheme dependence of the results turns out to be very large; hence, 
one can neither obtain accurate prediction of the potential in this distance region. For instance, 
the QCD potential bends downwards at large r as compared to the Coulomb potential if the 
l^-scheme running coupling constant is used, whereas the potential bends upwards at large r if 
the F-scheme running coupling constant is used [17]. (See e.g. Fig. 4 of [18].) It was later pointed 
out that the large uncertainty of the perturbative QCD prediction can be understood as caused 
by the C(Aq C d) infrared (IR) renormalon contained in Vqcd(^) [19]. 

Empirically it has been known that phenomenological potentials and lattice computations of 
VqcdM are both approximated well by the sum of a Coulomb potential and a linear potential in 
the above range 0.5 GeV" 1 < r < 5 GeV" 1 [20]. The linear behavior of VqcdM at large distances 
r 3> Aq^ d , verified numerically by lattice simulations, is consistent with the quark confinement 
picture. For this reason, and given the very poor predictability of perturbative QCD, it was often 
said that, while the "Coulomb" part of Vqcd(^) (with logarithmic corrections at short-distances) 
is contained in the perturbative QCD prediction, the linear part is purely non-perturbative and 
absent in the perturbative QCD prediction (even at r < Aq^ d ), and that the linear potential 
needs to be added to the perturbative prediciton to obtain the full QCD potential. Nevertheless, 

*Recently 2-loop correction to the octet QCD potential has been computed [16]. 
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to the best of our knowledge, there was no firm theoretical basis for this argument. 

Since the discovery [21, 22, 23] of the cancellation of C(Aq CD ) renormalons in the total energy 
of a static quark- ant iquark pair E tot (r) = VQCD( r ) + 2m polc , convergence of the perturbative series 
for E tot (r) improved drastically and much more accurate perturbative predictions for the potential 
shape became available. It was understood that a large uncertainty originating from the (9(Aqcd) 
renormalon in Vqcd(^) can be absorbed into twice of the quark pole mass 2m po i e - Once this 
is achieved, perturbative uncertainty of E tot (r) is estimated to be (9(AQ CD r 2 ) at r < Aq^ d [19], 
based on the renormalon dominance hypothesis. 

On the other hand, OPE of Vqcd(^) for r <C Aqq D was developed [10, 24] within an effective 
field theory "potential non-relativistic QCD" (pNRQCD) [25]. In this framework, Vqcd(?") is 
expanded in r (multipole expansion). At each order of this expansion, short-distance contributions 
are factorized into Wilson coefficients (perturbatively computable) and long-distance contributions 
into matrix elements of operators (non-perturbative quantities). The leading non-perturbative 
contribution to the potential is contained in the 0(r 2 ) term of the multipole expansion. 

Subsequently, several studies [18, 9, 26, 27] showed that perturbative predictions for Vqcd(^) 
agree well with phonomeno logical potentials and lattice calculations of Vqcd(^), once the O(Aqcd) 
renormalon contained in Vqcd(^) is cancelled. In particular, in the context of OPE, the leading 
Wilson coefficient was shown to be in agreement with lattice computations of Vqcd(^), after 
the subtraction of the (9(Aqcd) renormalon [26]. Ref. [28] showed that a Borel resummation of 
the perturbative series gives a potential shape which agrees with lattice results, if the O(Aqcd) 
renormalon is properly taken into account. In fact, these agreements hold within uncertainties of 
(9(Aq CD r 2 ) estimated from the residual renormalon. That is, a linear potential of C(AQ CD r) at 
r < Aq^ d was ruled out numerically in the differences between the perturbative predictions and 
phenomeno logical potentials/lattice results. These observations support the validity of renormalon 
dominance hypothesis. 

A crucial point is that, once the C(Aqcd) renormalon is cancelled and the perturbative pre- 
diction is made accurate, the perturbative potential becomes steeper than the Coulomb potential 
as r increases. This feature is understood, within perturbative QCD, as an effect of the running 
of the strong coupling constant [29, 18]. 

Soon after, it was shown analytically [30] that the perturbative QCD potential approaches a 
"Coulomb" +linear form at large orders, up to an (9(AQ CD r 2 ) uncertainty. (Here and hereafter, the 
"Coulomb" potential with quotes represents a Coulombic potential with logarithmic corrections 
at short distances.) Higher-order terms were estimated by the large-/9o approximation or by RG 
equation and a scale-fixing prescription based on renormalon dominance hypothesis was used. 
The "Coulomb" +linear potential can be computed systematically via RG; up to NNLL, it shows 
a convergence towards lattice computations of Vqcd(^)- Furthermore, the "Coulomb" +linear 
potential was shown to coincide with the leading Wilson coefficient in the framework of OPE, up 
to an 0(r 2 ) difference [31]. 

In this paper, we perform a precise and solid analysis, on the basis of our previous works 
[30, 31]. This work extends our previous works in the following respects: 

• We incorporate a degree of freedom for varying renormalization scale into the analysis of 
[30]. In this way, the "Coulomb" -l-linear potential is identified with the scale-independent 
part of the prediction. Details of the derivation and formulas not delivered so far are also 
presented. 

• We promote the "Coulomb" -l-linear potential to the leading Wilson coefficient in the frame- 
work of OPE, taking advantage of the result of [31]. We study properties of the Wilson 
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coefficient and the corresponding non-perturbative correction 5E\js(r). 
In addition, we present the following analysis: 

• We determine the non-perturbative correction 5Eus(r) using perturbative computations of 
the Wilson coefficients and recent lattice data. 

• As a byproduct, we determine the relation between lattice scale (Sommer scale) and A^g. 
This provides a new method to determine this relation. 

In this analysis, we assemble all the developments of perturbative computations and of OPE up 
to date. 

Organization of the paper is as follows. Sec. 2 is devoted to a review: we review the cur- 
rent status of the perturbative QCD computations of Vqcd(?") (Sec. 2.1), convergence property of 
Etot{ r ) U P to 0{a%) (Sec. 2.2), large-order behavior of the perturbative series based on renormalon 
argument (Sec. 2.3), and the predictions of OPE for Vqcd(^) (Sec. 2.4). In Sec. 3, we analyze 
the large-order behavior of the perturbative prediction of Vqcd(^) analytically: After explaining 
the strategy in Sec. 3.1, we present the results when the higher-order terms are estimated by the 
large-/3o approximation and by RG in Sees. 3.2 and 3.3, respectively. Details of the derivation 
are given through Sees. 3.4 and 3.5. (The readers may as well skip these details in the first read- 
ing.) Sec. 4 defines two types of renormalization schemes for the leading Wilson coefficient in the 
context of OPE (Sees. 4.1 and 4.2) and discusses properties of the Wilson coefficient and of the 
corresponding non-perturbative contributions (Sec. 4.3). In Sec. 5 we compare the perturbative 
computations of the Wilson coefficient with lattice compuations of Vqcd(?")- We first check consis- 
tency of theoretical predictions based on OPE (Sec. 5.1). Then we determine the non-perturbative 
contribution in each scheme as well as the relation between lattice scale and A^jg (Sees. 5.2 and 
5.3). Summary and conclusions are given in Sec. 6. 

App. A collects the formulas necessary for the computation of the perturbative series of the 
QCD potential. In App. B, we give a derivation of the one-parameter integral representation 
of [tty T (5')]oo- In App. C, we present the analytic formula for the linear potential up to NNLL. 
Methods for numerical evaluation of the Wilson coefficient are given in App. D. 

2 Perturbation Series and OPE of Vqcd(^) (Review) 

2.1 Definitions and conventions 

Throughout this paper, color factors of QCD are denoted as 



where Nc is the number of color, Cf is the second Casimir operator of the fundamental repre- 
sentation, Ca is the second Casimir operator of the adjoint representation, and Tp is the trace 
normalization of the fundamental representaion of the color SU (3) group. Furthermore, we denote 
the number of light quark flavors by n\. We assume that all light quarks are massless (except in 
Sec. 2.2). 

The static QCD potential is defined from an expectation value of the Wilson loop as 



4 1 

Cf — g, C A = N c = 3, Tp = -, 



(1) 



^ocd(^) = — lim — log 



(0\TrPexp[ig s f v dxi i A li (xj\ |0> 



(2) 



(0 | Tr 1 | ) 



3 



where P is a rectangular loop of spatial extent r and time extent T. The second line defines the 
^-scheme coupling contant, cav(q), in momentum space. In dimensional regularization, there are 
one temporal dimension and d = D — 1 = 3 — 2e spatial dimensions. 

In perturbative QCD, ay(q) is calculable in series expansion of the strong coupling constant. 
We denote the perturbative evaluation of ay(q) as 

a v T (q) = a 5 (/i)^P n (logWg))(^) (4) 

n=0 

= rtEwff) (5) 

n=0 

Here, as{fi) denotes the strong coupling constant renormalized at the renormalization scale fi, de- 
fined in the modified minimal subtraction (MS) scheme; P n (£) denotes an n-th-degree polynomial 
of £. In the second equality, we set \i = q using /i-independence of a v T (q). Eq. (5) is reduced to 
eq. (4), if we insert the series expansion of as(q) in terms of as(/z). This expansion is determined 
by the RG equation 

q 2 ^- 2 as(q)=(3(a s (q)) = -a s (q) £ /^(^f^ ( 6 ) 

n=— 1 

where f3 n represents the (n + l)-loop coefficient of the beta function.* Eqs. (4) (5) show that, at 
each order of the expansion of a v T (q) in ag(/x), the only part of the polynomial P n (\og(fj,/q)) that 
is not determined by the RG equation is P n (0). 

It is known [1] that P n (0) for n > 3 contain IR divergences. Namely, the perturbative QCD 
potential is IR divergent and not well-defined at and beyond 0(a s ). There are two ways to 
deal with this problem. One way is to use OPE, in which the QCD potential is factorized into 
Wilson coefficients and matrix elements. The Wilson coefficients include only ultraviolet (UV) 
contributions, hence they are computable in perturbative expansion in 0:5 free from IR divergences. 
IR contributions are contained in the matrix elements which are non-perturbative quantities. 
Another way is to expand the QCD potential as a double series in as and log as- This is achieved 
by resummation of certain class of diagrams (Fig. 1) as indicated by [1]. More systematically, 
this can be achieved within pNRQCD framework [10, 11, 24]. We will use both methods for 
regularization of IR divergences through Sees. 3-5. * 

Let us explain our terminology for the order counting. When we state u ay T (q) up to 0(ag)" 
we mean that we truncate the series on the right-hand-side of eq. (4) and take the sum for 
< n < N — 1. We also improve the perturbation series using the RG evolution of the MS 
coupling. * By a v T (q) up to LL, NLL, NNLL and NNNLL, we mean that we define ay T (q) by 

*In dimensional regularization and MS scheme, /3_i 7^ when the space-time dimension is different from 4; see 
App. A, eq. (123). 

tin our analysis in Sec. 3, regularization of IR divergences is rather a conceptual matter; there, our practical 
analysis concerns only up to the orders where IR finite terms are involved. On the other hand, in Sees. 4-5, we 
include 0{a A s ) term in our analysis, hence the regularization becomes practically relevant. 

tit is known that, up to NNLL, the RG-improved MS running coupling is more convergent than the RG-improved 
running coupling in the F-scheme or F-scheme, hence the RG-improvement in the MS-scheme leads to a more stable 
prediction of the potential shape; see [32] and Sec. 4 of [18]. For this reason, we adopt the RG-improvement in the 
MS-scheme in this paper. 
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Figure 1: Class of diagrams contributing to the QCD potential at O (a s log as). Dashed lines represent 
Coulomb gluons; curly line represents transverse gluon. 

eqs. (5) (6) and take the sums for < n < 1, 2, 3 and 4, respectively, in both equations (i.e. 1-, 2-, 
3- and 4-loop running coupling constants are used for as(q), respectively). This procedure resums 
logarithms of the forms as(fi)[as(f^) log(/i/g)] n , • • •, a,s(/z) 4 [a;s(/i) log(/z/g)] n , respectively. 

On the other hand, the IR divergences at 0(a s ) and beyond induce additional powers of 
log(/i/q) in otv(q) at NNLL and beyond, which are not resummed by the evolution of ots(q) via 
eq. (6). Hence, at these orders, it is more consistent (with respect to naive power counting) to 
resum these IR logarithms (referred usually as ultrasoft logarithms) as well, although physical 
origins of the logarithms are quite different. The ultrasoft (US) logarithms at NNLL can be 
resummed by replacing the U-scheme coupling constant as [12] 



a 



PT 

V 



(q)^a^(q) + ^a s (qf 



log 



a s (nf) 



(7) 



where \x$ denotes the factorization scale. We will examine the resummation of US logs separately. 

For n < 2, we define a n = P n (0). For n > 3, we include US logs into a n in addition. Explicit 
expressions for P n (£), a n , f3 n up to n = 3 (except for the unknown part of a 3 ) are listed in App. A. 
Furthermore, for convenience, we will denote 



S = ^1(31 



(8) 



in the following. Other formulas, useful for evaluation of oty(q), are collected in App. A as well. 



2.2 Convergence and scale-dependences of E tot (r) up to O(cx^) 

Let us demonstrate the improvement of accuracy of the perturbative prediction for the total 
energy E tot (r) = 2m polc + Vq C d(^) up to 0(a 3 s ), when the cancellation of (9(A QCD ) renormalons 
is incorporated. This is achieved (even without any knowledge of renormalons) if one re-expresses 
the quark pole mass m po i e by the MS mass in series expansion in as{fi). Presently perturbation 
series of VqcdM [4, 5] and m pole [33] are both known up to 0(a s ). 

As an example, we take the bottomonium case:§ We choose the MS mass of the 6-quark, 



renormalized at the 6-quark MS mass, as m,b 
flavors of light quarks are included with m u 
formula for E^ t (r) in [9].) In Fig. 2, we fix r : 



MS (m MS) 



2.5 GeV" 



= 4.190 GeV; in internal loops, four 
= and m c = 1.243 GeV. (See the 
0.5 fm (midst of the distance range 



of our interest) and examine the renormalization scale (//) dependence of E tot (r). We see that 
E tot (r) is much less scale dependent when we use the MS mass (after cancellation of renormalons) 
than when we use the pole mass (before cancellation of renormalons). This shows clearly that the 
perturbative prediction of E tot (r) is much more stable in the former scheme. 

§E^ t (r) = 277Jb,poic + VqcdW has been applied to computations of the bottomonium spectrum [34]. 
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Figure 2: Scale dependences of E^ t (r) up to C(a|) at r = 2.5 GeV 1 0.5 fm, in the pole-mass and 
MS-mass schemes. A horizontal line at 8 GeV is shown for a guide. 

We also compare the convergence behaviors of the perturbative series of E tot (r) for the same 
r and when fj, is fixed to the minimal-sensitivity scale [35] (the scale at which E tot becomes least 
sensitive to variation of /i) in the MS-mass scheme. At r = 2.5 GeV -1 , the minimal-sensitivity 
scale is \i = 0.90 GeV. Convergence of the perturbation series turns out to be close to optimal for 
this scale choice:^ 



The four numbers represent the 0(afg), 0{a\), 0(a 2 s ) and 0{a i s ) terms of the series expansion 
in each scheme. The O(a s ) terms represent the twice of the pole mass and of the MS mass, 
respectively. As can be seen, if we use the pole mass, the series is not converging beyond 0(ag), 
whereas in the MS-mass scheme, the series is converging. One may further verify that, when the 
series is converging (MS-mass scheme), //-dependence of E tot (r) decreases as we include more terms 
of the perturbative series, whereas when the series is diverging (pole- mass scheme), /i-dependence 
does not decrease with increasing order. (See e.g. [36].) 

We observe qualitatively the same features at different r and for different number of light 
quark flavors ni, or even if we change values of the masses m^, m c . Generally, at smaller r, E tot (r) 
becomes less /i-dependent and more convergent, due to the asymptotic freedom of QCD [9]. 

The stability against scale variation and convergence of the perturbative series are closely 
connected with each other. Formally, scale dependence vanishes at all order of perturbation series. 
This means that, for a truncated perturbative series up to 0(aig), scale dependence is of 0{a^ +1 ). 
Hence, the scale dependence decreases for larger N as long as the series is converging. Thus, the 
truncated perturbative series is expected to become less //-dependent with increasing order when 
the series is converging. It also follows that the series is expected to be most convergent when 

'in the pole-mass scheme, there exists no minimal-sensitivity scale within a wide range of /i, and the convergence 
behavior of the series is qualitatively similar to eq. (9) within this range. 




10.408 - 0.275 - 0.362 - 0.784 GeV (Pole-mass scheme) 
8.380 + 1.560 - 0.116 - 0.022 GeV (MS-mass scheme). 



(9) 
(10) 
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[i is close to the minimal-sensitivity scale. This observation, supported by the above numerical 
verification up to forms a basis of our analysis in Sec. 3. 

As already mentioned in the Introduction, once E tot (r) is expressed in terms of the MS mass 
and an accurate prediction is obtained, it agrees well with phenomenological potentials and lattice 
computations of the QCD potential in the range of r of our interest. As more terms of the series 
expansion are included, E tot (r) becomes steeper in this range. This behavior originates from an 
increase of the interquark force due to the running of the strong coupling constant [18]." E to t(r) up 
to a finite order in perturbative expansion has a functional form 1/rx (Polynomial of logr), apart 
from an r-independent constant; cf. App. A, eq. (126). On the other hand, we see a tendency that, 
as we increase the order, E tot (r) approaches phenomenological potentials/lattice results, which are 
typically represented by a Coulomb+linear potential. This observation motivates us to examine 
the perturbative prediction for E tot (r) at large orders, which will be given in Sec. 3. For that 
analysis, we need to know large-order behaviors of the perturbative series of E tot (r). 

2.3 Large-order behaviors and IR renormalons 

The nature of the perturbative series of VqcdM and E tot (r) at large orders, including their 
uncertainties, can be understood within the argument based on renormalons. The argument gives 
certain estimates of higher-order terms, and empirically it gives good estimates even at relatively 
low orders of perturbative series. 

Before starting any argument on large-order behaviors, one may be perplexed because the 
perturbative expansion of Vqcd(^) contains IR divergences beyond 0{a\). For definiteness, let 
us assume (conceptually) that we regularize the IR divergences by expanding VqcdM in double 
series in a s and log a s ; then we identify O(a^) term with the sum of 0(ag\og k a s ) terms for all 
k* 

Let us denote the 0{a^ +1 ) term as V^! D (r). According to the renormalon argument, the 
leading behavior of Vq^ D (r) at large orders n ^> 1 is given by 

KW(r)~ const, x („) 

up to a relative correction of 0(l/n) [38]. It follows that (Vg^^r)! becomes minimal at order 
n ~ N = 27r/((3 as(iJi)), while |VQ^ D (r)| scarcely changes in the range N — V^Ao <C n 
No + y/No- For n ^> iVo + V^o> the series diverges rapidly. (See Fig. 3, black squares.) Due to the 
divergence (the series is an asymptotic series), there is a limitation to the achievable accuracy of 
the perturbative prediction for Vqcd(^)- An uncertainty of the asymptotic series may be estimated 
by the size of the terms around the minimum, a/Aq x | ^qcd ( r ) I ' which gives an uncertainty of 
O(Aqcd) [19]. 

The perturbative series of E tot (r) in the pole- mass scheme is the same as that of Vqcd(^) 
except for the O(a s ) term. If we re-express E tot (r) in terms of the MS mass, the leading behavior 
of Vqq D (t) is cancelled against that of the perturbative series of 2m vo \ e ) Then the large-order 

"See [29, 37] for a more microscopic explanation of this feature. 

* There exists evidence that renormalon dominance may be valid in such an expansion [36]. 

tin order to realize the cancellation of the leading behavior of the perturbative series at each order of the 
expansion, one needs to expand Vqct>{t) and m po i in the same coupling constant otsip). This is somewhat 
involved technically, since usually VqcdM an d m po \ c are expressed in terms of different coupling constants; see 
[18, 9, 27]. 
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Figure 3: Diagram showing the n-dependence of |Vqq D ( 



(or 1 42 (r) 



squares] and that of \E, 



in). 



in the pole-mass scheme) [black 



in the MS-mass scheme [red squares], based on renormalon estimates. 



behavior of E tot {r) becomes 

E${r) ~ const, x r 2 n\ (^^Y (12) 

|i?tot( r )l becomes minimal at n ~ Aq = 67r/(/3 tts(/i)) an d its size scarcely changes for Aq — 
y/Ni C C iVj + y/Ni- As compared to E tot (r) in the pole-mass scheme, the series converges 
faster and up to a larger order, but beyond order a^ 1 again the series diverges. (See Fig. 3, red 
squares.) An uncertainty of the perturbative prediction for E tot (r) can be estimated similarly as 
v^xl^^l-^A^r 2 ) [19]. 

We note that each term of the perturbation series {V^ D , E^) is dependent on the scale /i. 
Hence, its large-order behavior, including the order at which its size becomes minimal [Nq,Ni oc 
l/ats(p)], is a l so dependent on \x. The estimated uncertainty (^/N^x \V$$(r)\, y/N^x \E^\r)\), 
however, is independent of /i. 

These estimates of large-order behaviors, according to renormalons, follow primarily from 
analyses of IR sensitivities of certain classes of Feynman diagrams; then the estimates are improved 
and reinforced via consistency with RG equation [38]. The C(Aq CD ) IR renormalon, corresponding 
to the perturbative series 

c r n en (k) ~ const, x n\ ^g-^T j" n kS , (13) 
originates typically from an integral of the form 

/ dqq 2k - 1 a s ( q ) = J2cn n (k), (14) 

where fif ^> Aqcd is a UV cutoff. Nevertheless, in general contributions originate also from more 
complicated loop integrals. 
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Q < Hf 




g s r- E a 



Figure 4: Leading contribution of US gluon to SEjjs(r) in pNRQCD. 



2.4 OPE of V QC n(r) 

A most solid way to separate perturbative and non-perturbative contributions to the QCD po- 
tential is to use OPE. OPE of the QCD potential was developed [10, 24] within pNRQCD [25], 
which is an effective field theory (EFT) tailored to describe dynamics of ultrasoft gluons coupled 
to a quark- ant iquark (QQ) system, when the distance r between Q and Q is small, and when 
the motions of Q and Q are non-relativistic. (In the case of the QCD potential, they are static.) 
Within this EFT, the QCD potential is expanded in r (multipole expansion), when the following 
hierarchy of scales exists: 

A QC D</i/<^ (15) 

Here, /if denotes the factorization scale. Non-perturbative contributions to the QCD potential 
are factorized into matrix elements of operators, while short- distance contributions are factorized 
into potentials, which are in fact Wilson coefficients. Conceptually, physics from IR region q < /if 
is contained in the former, while physics from UV region q > fif is contained in the latter. 
Explicitly, the QCD potential is given by [24] 

V Q cD(r) = V s (r) + 5E lJS (r), (16) 

5E vs (r) = -igl^r dte- iAv ^U0\r-E a (t)^ di (t,0) ab r-E b (0)\0) + 0(r 3 ). (17) 
Jo 

The leading short-distance contribution to Vqcd(^) is given by the singlet potential Vs(r). It is a 
Wilson coefficient, which represents the potential between the static QQ pair in the color singlet 
state. The leading long-distance contribution is contained in the matrix element in eq. (17). It is 
0(r 2 ) in multipole expansion. AV(r) = Vo(r) — Vs{r) denotes the difference between the octet 
and singlet potentials; E a denotes the color electric field at the center of gravity of the QQ system. 
See [24] for details. 

Intuitively we may understand why the leading non-perturbative matrix element is 0(r 2 ) as 
follows. As well known, the leading interaction (in expansion in r) between soft gluons and a 
color-singlet QQ state of size r is given by the dipole interaction r ■ E a . It turns the color singlet 
QQ state into a color octet QQ state by emission of soft gluon(s). To return to the color singlet 
QQ state, the color octet state needs to reabsorb the soft gluon(s), which requires an additional 
dipole interaction. Thus, the leading contribution of soft gluons to the total energy is 0(r 2 ). See 
Fig. 4. 

Although 5E\js(r) is 0(r 2 ) in terms of the expansion of operators, it has an additional depen- 
dence on r through the Wilson coefficient AV(r). After all, we would like to know how SE\js(r) 
depends on r in the region of our interest. The leading power of r can be determined in some 
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cases. Since, however, the argument depends on the renormalization of the singlet potential within 
pNRQCD, let us discuss this issue first. 

The Wilson coefficient Vs(r) can be computed in perturbative expansion in as by matching 
pNRQCD to QCD. It turns out that Vg(r) thus computed coincides with the perturbative expan- 
sion of Vqcd(^) (in dimensional regularization); in particular, this means that Vs(r) includes IR 
divergences beyond 0(a s ). This result follows from a simple argument: Formally, SE\js(r) can be 
computed also in series expansion in as- This expansion, in dimensional regularization, vanishes 
to all orders, since all diagrams are given by scaleless integrals.* 

On the other hand, SE\js(r) is expected to be non-zero beyond naive perturbation theory. For 
instance, this can be verified by computing 5E\js(r) in pNRQCD when as (1/r) <C 1. According 
to the concept of the EFT, Vs(r) and AV(r) should be expanded in as only after all loop in- 
tegrations are carried out. Since this theory is assumed to correctly describe physics at energy 
scales much below 1/r, AV(r) (^C 1/r) should be kept in the denominator of the propagator 
[E — Al^r)] -1 .* Thus, if we expand all factors except AV(r) in 0:5 in eq. (17), 5Eus(r) becomes 
non-zero since AV(r) acts as an IR regulator. (One may expand AV(r) in as only after all the 
integrations are performed. Then log as appears, in contrast to the formal expansion in as, where 
everything is expanded before integrations.) In this case, 5E\j$(r) contains UV divergences. In 
dimensional regularization (D = 4 — 2e), they are given as poles in e, which exactly cancel the 
poles corresponding to the IR divergences in Vs(r). Consequently, in the sum eq. (16), Vqcd(^) 
becomes finite as e — > 0. 

These divergences in Vs(r) and 5E\js(r), respectively, can be regarded as artefacts of dimen- 
sional regularization, where the integral regions of virtual momenta extend from to 00. If we 
introduce a hard cutoff to each momentum integration, corresponding to the factorization scale 
[if, Vs(r) (q > fif) and 5Eu S (r) (q < fif), respectively, would become finite and dependent on fif. 
This observation calls for renormalization of Vs(r) and 5E\js(r) within pNRQCD also in dimen- 
sional regularization. For example, Vs(r) can be made finite by multiplicative renormalization, 
i.e. by adding a counter term (Zs — l)Vs(r). 

With respect to the spirit of factorization in OPE, it is natural to subtract IR renormalons 
from Vs{r) in a similar manner. In [14, 26], this was advocated and in practice subtraction of 
(only) the O(Aqcd) renormalon was carried out explicitly. The known IR renormalons of the bare 
Vsir) are contained in the integral [38]* 



[Note that the perturbative expansion of the bare Vs(r) coincides with that of VqcdM-] As for 
the 0(Aq CD r 2 ) renormalon, it was shown that the IR renormalon contained in the bare Vs(r) and 
the UV renormalon contained in the bare 5E\js(r) cancel in dimensional regularization [31]. In 
a hard cutoff renormalization scheme, contributions of gluons to 5Eu S (r) close to the UV cutoff 
region q ~ /if can be analyzed using perturbative expansion in as within pNRQCD, due to 
the hierarchy (15). It has exactly the structure suitable to absorb the (9(AQ CD r 2 ) renormalon 

*We neglect the masses of quarks in internal loops. 

TThis situation is similar to the case, where one should not expand the electron propagator by the electron mass 
if one wants to describe the physics of collinear photon emission in the region Ed <C m e . 

^Here, we neglect the contributions of the instanton-anti-instanton-induced singularities [38] on the positive real 
axis in the Borel plane. These contributions are known to be rather small. 




(18) 
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contained in eq. (18). Namely, in a hard cutoff scheme, the 0(AQ CD r 2 ) renormalon is subtracted 
from V s (r) and absorbed into 5E vs (r). The /independences that enter as a consequence cancel 
between the renormalized Vs(r) and 5E\js(r) [24]. Hence, everything holds in parallel with the 
case of IR divergences discussed above. Therefore, it is appropriate to subtract from Vs(r) the IR 
renormalons, e.g. in the form of eq. (18), in addition to subtracting IR divergences, and to define 
a renormalized singlet potential. (We will give explicit renormalization prescriptions in Sec. 4.) 

More generally, it is known that, in a wide class of physical observables (whenever OPE is avail- 
able), IR renormalons in perturbation series are deeply connected with OPE of the corresponding 
physical observables. As we have seen, renormalon uncertainties have power-like behaviors in the 
ratio of a large scale and Aqcd [i n our case (r Aqcd)* 1 x Aqcd]- in OPE, non-perturbative con- 
tributions (matrix elements of operators) have the same power-like structures. Therefore, in an 
appropriate renormalization prescription, IR renormalons contained in perturbative series can be 
subtracted from Wilson coefficients and absorbed into matrix elements in OPE, thereby leaving 
Wilson coefficients free from IR renormalons. It means that (in principle) Wilson coefficients can 
be computed to arbitrary accuracy by perturbative expansion. At the same time, renormalon 
ambiguities are replaced by matrix elements of operators (condensates), the values of which can 
be determined by comparing to various experimental data or results of lattice simulations. 

Now we return to the discussion on the r-dependence of SEu S (r) when r <C Aq CD [24]. We 
assume that V^(r) and 5E\js(r) are renormalized in a hard cutoff scheme, according to the above 
discussion. One can derive the r-dependence of SE\js(r) clearly when AV(r) « Caocs/t ^ A 4 / (^ 
Aqcd)- Since, in this case, the exponential factor in eq. (17) is rapidly oscillating, we can expand 
the matrix element in t. Then the matrix element reduces to a local gluon condensate, and from 
purely dimensional analysis, 5Eu S (r) becomes 0(fj,jr 3 ).§ The condition AV(r) 3> /!/ 3> Aq CD is 
satisfied at sufficiently short distances. 

Another case, in which r-dependence of 5E\js(r) is known, is when fif ^> AV(r) is satisfied, 
in addition to the hierarchy (15). This condition is expected to hold at r Aq CD but not for too 
small r. Under this condition, 5E\js(r) is dominated by contributions of gluons from the region 
Aqcd, AV(r) C g < /if, which can be computed in perturbative expansion in as- This leads to 
SE vs (r) ~ (?i,,}r-). 

Let us discuss the case where ///is reduced and taken close to Aqcd- This case violates the 
conventional hierarchy condition (15). If AV(r) ^> Aqcd, the matrix element can still be reduced 
to the local gluon condensate, and 5Eu S (r) ~ (9(AQ CD r 3 ). On the other hand, if AV(r) ~ Aqcd, 
there is no way to predict the r-depndence of 5Eu S (r) in a model-independent way. If AV(r) <C 
Aqcd, we can expand the exponential factor in AV(r) in eq. (17) and find 5E\js(r) ~ (9(AQ CD r 2 ).^ 

In the distance range of our interest, 0.1 fm < r < 1 fm, the relation between Aqcd and 
AV(r) is not very clear. A rough estimate shows that, at small r within this range (perhaps 
r < 0.3 fm), AV(r) 3> Aqcd, whereas at larger r (perhaps r > 0.3 fm), AV(r) ~ Aq C d- However, 
of course, this depends on a precise definition of Aqcd and accurate knowledge of AV(r). It is 
quite probable that there exists no [xj within the above range of r such that AV(r) ^> fif ^> Aqcd 
can be satisfied. Therefore, if we choose iif ^> Aqcd, we would expect 5Eus(r) ~ (9(/i 3 r 2 ) in 
the entire range 0.1 fm < r < 1 fm. On the other hand, if we choose \Xf ~ Aqcd, we conjecture 
that at small distances (perhaps 0.1 fm < r < 0.3 fm), 5E\j S (r) ~ (9(AQ CD r 3 ), whereas at larger 
distances, we cannot predict the r-dependence of 5E\js(r) in a model-independent way. 

§Note that we may ignore Aqcd in comparison to [if, since [if ^> Aqcd- An alternative derivation is to compute 
contributions of gluons from the region Aqcd q Hf using perturbative expansion in 0:5. 

^In this paper, we do not consider the possibility AV(r) <C Aqcd henceforth, since such large r seem to lie 
beyond the applicable range of our analysis. 
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To end this subsection, let us discuss what is indicated by OPE of the QCD potential as given 
above. Suppose we consider an expansion of Vqcd(?") at r < Aq^ d (in the distance range of our 
interest): 

VqcdM « ^ + Co + Cl r + c 2 r 2 + ---. (19) 

This is (at best) only a qualitative argument, since we know that there are logarithmic corrections 
to the Coulomb potential at short-distances, and for this reason, VqcdM cannot be expanded 
in Laurent series. Nevertheless, empirically the above expansion is a good one, since many phe- 
nomenological potentials have been successfully determined, by fitting them to the experimental 
data of heavy quarkonium spectra, assuming Coulomb+linear forms. So, suppose that one may 
decompose Vqcd(?") as above qualitatively. Then, since the non-perturbative contribution 5E\js(r) 
is expected to be 0(r 2 ) (assuming /// AV, Aqcd), the c 2 r 2 term (and beyond) would come from 
both Vs(r) and SEjjs(r), and their relative contributions change as we vary the factorization scale 
/if. On the other hand, the Coulomb, constant, and linear terms, c_i/r + c + Cir, should orig- 
inate only from the perturbative prediction of Vg(r), that is, from the perturbative prediction 
°f Vqcd(^)- (The constant term becomes predictable perturbatively only when the pole masses 
are added to Vqcd(?") and rewritten in terms of a short- distance mass such as the MS-mass.) In 
particular, they should be predictable independently oi/if. 

3 Perturbative QCD Potential at Large Orders 

In this section, we present an analysis of the QCD potential at large orders of perturbative ex- 
pansion. We separate the perturbative prediction of the QCD potential at large orders into a 
scale-independent (prescription-independent) part and scale-dependent (prescription-dependent) 
part, when higher order terms are estimated via large-/5o approximation or via RG, and when the 
renormalization scale \i is varied around the minimal-sensitivity scale. 

3.1 Strategy and general assumptions of the analysis 

We consider the perturbative QCD potential up to 0(ag): 

j^s -r K t (?)]at. (20) 

Here and hereafter, [X]n denotes the series expansion of X in cxs(fj) truncated at 0(as(fj) N ). We 
examine V N (r) for N >> 1. For this analysis, we need (a) an estimate for the all order terms of 
Vjv(r), and (b) a scale-fixing prescription. 

In the following subsections, we estimate the higher-order terms of V)v(r) using large-/3 ap- 
proximation (Sec. 3.2) and using RG (Sec. 3.3). There is a caveat: The former estimate does 
not contain IR divergences at all, and in the latter estimate, IR divergences appear only beyond 
NNLL; hence, in most of our argument, we will discard IR divergences. Since the true higher-order 
terms contain IR divergences beyond 0(a|), we have to clarify what we mean by our estimates 
of higher-order terms. Conceptually, we assume that we have removed ambiguities related to IR 
divergences, while keeping IR renormalons in the perturbative expansion of the potential. This 
seems to be possible, since, up to our current best knowledge, IR divergences [1] and IR renor- 
malons [19] contained in the perturbative QCD potential stem from quite different physical origins. 
As an explicit example to realize such a situation, we may assume that we analyze the singlet 
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potential Vg(r) instead of VqcdO"), after subtracting IR divergences (but not IR renormalons) 
via renormalization. Alternatively, we may assume that we have regularized IR divergences by 
expanding Vqcd(?") m double series in as and log as- 
Let us explain our scale-fixing prescription (b). Since within our estimates the perturbative 
series turns out to be an asymptotic series, there exists a certain arbitrariness in making a pre- 
diction from large-order analysis of the series. We will give a prediction by choosing a reasonable 
scale \i for each given N and then taking the limit iV — > oo. (Later we will justify our prescription 
by comparing the prediction with that in OPE.) Perturbative QCD in itself does not provide any 
scale-fixing procedure. In practice, whenever a perturbative expansion up to some finite order is 
given, one chooses a reasonable (range of) scale /i, as we have seen in Sec. 2.2. We would like 
to fix the scale in a similar manner in our large-order analysis. According to the argument given 
in Sees. 2.2 and 2.3, if we choose a scale \x such that a s (/i) = 6n/(p N) is satisfied, around this 
scale, Viv(r) (after cancelling the leading-order renomalon) would become least /x-dependent and 
the perturbative series would become most convergent; cf. Fig. 3. In view of this property, we fix 
fj, such that* 

JV = s£50 { = * t - (21) 

£ = 1 corresponds to an optimal choice; by varying the parameter £, we may change the scale 
ji for a given N . Then we consider V/v(r) for N ^> 1 while keeping Aj^g finite. (Here, we relate 
our scale-fixing prescription to that of principle of minimal sensitivity [35] only weakly, as argued 
above. A close examination of the relation can be found in [39].) 

An alternative way to regard this prescription is as follows. Suppose we know the perturbative 
expansion of Vqcd(?") up to all orders, according to a certain estimate. When the expansion is 
asymptotic for any choice of cxs(fj), we cannot sum all the terms. Instead, following a standard 
prescription to deal with asymptotic series, we may truncate the series around the order where 
the term is close to minimal. This gives the truncated series Vjv(r) with N given by the relation 
eq- (21). 

The motivation for considering the large N limit is that it corresponds to the limits where 
the perturbative expansion becomes well-behaved (small expansion parameter) and where the 
estimate of Vqq^t) by renormalon contribution becomes a better approximation around n ~ N. 
Note that large N corresponds to small as(n) and large \x due to the above relation. 

Let us further comment on some details concerning the relation (21). (i) The relation (21) 
follows from the asymptotic form, eq. (12), of the series independently of its overall coefficient. 
Although the overall coefficient is not known exactly,' other parts of eq. (12) or (13) are considered 
to be solid, based on consistency with RG equation. Hence, the relation (21) is based on a 
solid part of the renormalon estimate, (ii) The scale \x fixed by the relation (21) is independent 
of r. Usually it is considered that a natural choice of the scale is related to a physical scale, 
typically /i ~ 1/r, at low orders of perturbative expansion. Moreover, the minimal-sensitivity 
scales corresponding to low orders of perturbative expansion, as in the cases of Sec. 2.2, are 
known to be strongly dependent on r [18, 9]. This is, however, not expected to be the case at 
large orders. It is because, in eq. (18), contributions from q < 1/r are dominant on the left-hand- 
side at low orders, whereas at large orders, the term proportional to — q 2 r 2 /6 dominates on the 
right-hand-side of eq. (18), hence, r 2 factors out as an overall coefficient; cf. eq. (12). (iii) Based 
on the argument in Sec. 2.3, we may consider that an optimal choice of \x or £ corresponds to the 
range N x - JW X < N = < N t + y/N^ in the relation (21). Then, £ -> 1 as N -> oo. 

*Here, we generalize the prescription of [30] by introducing an additional parameter £, cf. [39]. 
'''See [40] for a method for systematically estimating the overall coefficient. 
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3.2 V/v(r) at large orders: large-/3 approximation 



The large-/3o approximation [41] is an empirically successful method for estimating higher-order 
corrections in perturbative QCD calculations; see e.g. [38, 33, 42, 43]. In general, the large-/3o 
approximation of a physical quantity, at a given order of perturbative expansion in as, is defined 
in the following way. We first compute the leading order contribution in an expansion in l/ni, 
which comes from so-called bubble chain diagrams. Then we transform this large n\ result by a 
simplistic replacement ni — > n\ — 33/2 = — (3/2)(3 . 

For the QCD potential, the large-/3 approximation corresponds to setting a n = (5j3 /3) n in 
eq. (5) and all (3 n = except (3q in eq. (6). Hence, it includes only the one- loop running of 
ots(q). In this subsection, with these estimates of the all-order terms of ay T (q), we examine 
Vjv(r), defined above, for AT ^> 1. The reasons for examining the large-/3 approximation are as 
follows. First, because this approximation leads to the renormalon dominance picture; in fact, 
the renormalon dominance picture has often been discussed in this approximation. Secondly, 
the running of the strong coupling constant makes the potential steeper at large distances as 
compared to the Coulomb potential; hence, we would like to see if the potential can be written in 
a "Coulomb" +linear form when only the one-loop running is incorporated as a simplest case. 
We define A 



e 5/ 6 A Woo P where 

MS ' 



A l-loop 

A MS =^ X P 



2tt 



In the following, we assume 



A" 



( N \ 7-1 ( N 

expl-— l<r<A exp( — 



(22) 



(23) 



oo. Note that, as A" — > oo, 



3^ r \H> 

when we consider the double limits r^0,Af^ooorr^oo,A r 
the lower bound (A _1 e _iV ^ 3 ^) and the upper bound (A~ 1 e 7V ^ 3 ^) of r go to and oo, respectively. 

First we present the result and discuss some properties when £ = 1, which corresponds to an 
optimal choice of scale p. (Derivation will be given in Sec. 3.4.) 



Result for £ = 1 

Vn(t) for £ = 1 and A^ ^> 1 within the large-/3 approximation can be decomposed into four 
parts corresponding to {r _1 , r°, r 1 , r 2 } terms (with logarithmic corrections in the r _1 and r 2 
terms) : 



4C P 



A v(Ar, N), 



(24) 



v (p, N) = v c (p) + B(N) + Cp + D(p, N) + (terms that vanish as A^ -> oo). (25) 
(i) "Coulomb" part: 

tt/2 



7T 1 [°° 

vc{p) = 1 — / dxe x arctan 

P PJo 



\og(p/x) 



(26) 



where arctan x G [0, 7r). The asymptotic forms are given by 



vc(p) 



7T 



vc(p) ~ - -, 

p 



2plog(l/p)' 



7T 



p — > oo 



(27) 
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and both asymptotic forms are smoothly interpolated in the intermediate region. The short- 
distance behavior is consistent with the one-loop RG equation for the QCD potential, 
(ii) constant part*: 



B(N) = 



-t r 



dt- 



(l + lt 



N 



N 



log 2 



9 99 
8N + 6AN 2 ' 



(28) 



The first term (integral) diverges rapidly for N 
(iii) linear part: 




i + e>(i/iv)p 



(iv) quadratic part: 



=p 2 


- 1 

.12 


POO 


< 


L 


dx - 








1 

12 


log( 



log N + d{p) 

( ■* — 1 — x + \x 2 — \x 3 6(1 — x) 



(29) 



(30) 



\og(p/x) 



x> 



\0g\ P /x) + 7T 2 /4 



log(log 2 p + + log ^ + 7s 



(31) 



where 6>(x) is the unit step function and ^ E = 0.5772... is the Euler constant. The asymptotic 
forms of d(p) are given by 



1 r 9 
d(p) ~ [21oglog(l/p) +log- + 7s 



1 r 9 
~ - ^2 [ 2 log log P + log - + 7s 



p^O 



p — > oo 



(32) 



and in the intermediate region both asymptotic forms are smoothly interpolated. 

The "Coulomb" [i>c(p)], linear [Cp] and quadratic [D(p, N)] parts are shown in Fig. 5. The 
truncated potential v(p,N) is compared with the "Coulomb" +linear potential vc(p) + Cp after 
the constant B(N) is subtracted, for iV = 10, 30, 100, in Fig. 6.* In order to show how quickly 
v(p,N) approaches v c (p) + B(N) + Cp + D(p,N) as iV increases, we show their differences for 
several values of iV in Fig. 7. One sees that convergence is quite good at p < 1 (r < A -1 ) for 
N > 10. (For the purpose of separating different lines visibly, we plot potentials up to fairly large 
distances in this section. Nevertheless, we stress that, in most cases, our interests are in the region 
r < A-M) 

Although the constant part of vjf°\r)\g =1 diverges rapidly as N — > oo, the divergence can 
be absorbed into the quark masses in the computation of the total energy E tot (r) (or the heavy 
quarkonium spectrum). Therefore, in our analysis, we will not be concerned with the constant 
part of the potential but only with the r-dependent terms. 

*The 0(1/ N) and 0{l/N 2 ) terms in eq. (28) are irrelevant for N -> oo. We keep these terms in B(N) for 
convenience m examining Vjf°\r) at finite N; see Fig. 6 below. 

I'The integral can be expressed in terms of confluent hypergeometric function. 

*One can find formulas convenient for computing Vn{t) for a finite but large N in App. A. 

^Roughly speaking, one may regard A -1 ~ f fm. 
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Figure 6: Truncated potential after the constant term is subtracted, v (p, N) — B(N), (dashed) vs. p for 
N = 10, 30, 100 and £ = 1. "Coulomb" +linear potential, vc(p) + Cp, (solid black) is also plotted, which 
is hardly distinguishable from the N = 30 curve. 
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The quadratic part of Vjf°\r)\^ =1 diverges slowly as A 3 r 2 logiV ~ A 3 r 2 loglog(p/A). The 
dependence of Vjf (r) on A" is mild (after the constant part is subtracted); for instance, as shown 
in Fig. 6, the variation of v(p, N) — B(N) is small in the range r < A -1 as we vary A" from 10 to 
100; it corresponds to a variation of /i/A^ op from 30 to 3 x 10 14 . 



The "Coulomb" part and the linear part are finite as A" — > oo. In Fig. 6, we see that V 1 



(A>) ( 

N I 



is approximated fairly well by the sum of the "Coulomb" part and the linear part (up to an r- 
independent constant) in the region r < A -1 when we vary A" between 10 and 100. Moreover, as 
long as j2 log A" < 0(1), the difference between Vjf ' (r)L =1 and the "Coulomb" +linear potential 

remains at or below (9(A 3 r 2 ) in the entire range of r. Note that v(p, N) in this figure have the 
form of 1/r x (Polynomial of logr), and a priori it is not obvious at all that they approximate a 
"Coulomb" +linear potential. 

Results for £ ^ 1 

We vary £ in the scale-fixing prescription eq. (21) and decompose Vjf°\r) as in eqs. (24) and 
(25). As a salient feature, we obtain the same "Coulomb" +linear potential, vc(p) + Cp, as in the 
£ = 1 case. On the other hand, the constant B(N) and D(p,N) change. The latter no longer 
takes a quadratic form. Let us list how D(p,N) change with £. (See Fig. 8.) 

• 2/3 < £ < 1, 

D(p, N) is finite as A" — > oo: 



x e x — ( 1 — x + \x 2 



-37ri§/2 



\og(p/x) — iir/2 



(33) 



D(p, oo) = -p* x I dx v 2 J Im 

Jo 

Its asymptotic forms are given by 

Dip, oo) ^p 3 *- 1 — L- x r(-3£)sinf-7r£Y p^0 or p -> oo. (34) 
logp V2 / 

The asymptotic forms at p — > and p — >• oo have opposite signs. In the intermediate region 
D(p, oo) changes sign once. D(p, oo) for £ = 0.85, 0.9, 0.95 are plotted in Fig. 9. 
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Figure 8: v(p, N) for different values of £ and TV. (Dashed lines for £ = 0.9 and dot-dahed line for £ = 1.1.) 
For comparison, the "Coulomb" +linear potential vc(p) + Cp is also shown (solid line). Constants have 
been added to v(p, N) to make them coincide with vc(p) + Cp at p = 0.5. 
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• £>1 



Even powers of p, corresponding to IR renormalons, become more divergent as we increase 

D(p, N) = d 2 (N) p 2 + d 4 (N) p 4 + • • • + d v {N) p v + (finite term as N -> oo), (35) 

where z/ is the largest even integer satisfying v < 3£ — 1. di(N) diverges at least loga- 
rithmically (typically exponentially) as N — > oo. It diverges more rapidly for larger £ and 
smaller %. The asymptotic form of the finite (^-independent) term as p — > or p — > oo is 
p3?-i x Q Q g correc tion). 

• t < 2/3, 

_D(p, iV) becomes more dominant than the linear potential Cp at short-distances. We do not 
consider this possibility henceforth. (^ = 2/3 is marginal; the asymptotic form eq. (34) is 
valid at p — > but not at p — > oo.) 

Dependence of B(N) on £ is similar: It diverges more rapidly as iV — > oo for larger £, while it 
becomes finite when £ < 1/3. 

Thus, B(N) and D(p,N) are dependent on £, i.e. on the choice of scale via eq. (21); they are 
also divergent as iV — > oo for a sufficiently large £. Namely, B(N) and D(p,N) are dependent 
on the prescription we adopted to define our prediction. It is natural to consider the prescription 
dependence as indicating uncertainties of our prediction. In fact, B(N) and D(p, N) are associated, 
respectively, with the O(Aqcd) IR renormalon and (9(AQ CD r 2 ) IR renormalon (and beyond) in 
^qcd(^)- We have already seen that these renormalons induce uncertainties. On the other hand, 
the "Coulomb" +linear part [vc{p) + Cp] are independent of £ and N. Hence, vc(p) + Cp can be 
regarded as a genuine part of the prediction. In this regard, we remind the reader that there are 
no IR renormalons associated with the 1/r and r terms in the QCD potential [19]. 

One may associate the C(AQ CD r 2 ) renormalon with D(p,N) through following observations. 

(1) When £ = 1, the quadratic part of vjf°\r) diverges as A 3 r 2 log N. If the series expansion of 
m poic( m MS' a s) 01 ^qcd(^) is truncated at the order corresponding to the minimal term of the 
LO renormalon contribution, i.e. at order Nq = 2n/ (p as), [^ P oic]Af or f^QCD^)]^ diverges as 
AlogAo within the large-/3o approximation. We may compare AlogAf with the usual interpre- 
tation that m po ie and Vqcd(^) contain (9(Aqcd) perturbative uncertainties due to the LO renor- 
malons. (2) An argument similar to (1) applies for £ ^ 1. (3) As we will see in the next subsection, 
even if we estimate higher-order terms using RG equation and incorporate effects of the two-loop 
running and beyond, D(p, N) has a similar behavior to that in the large-/3o approximation. 

Let us further discuss questions concerning the strategy and results of the analysis given above. 

Naively, one would expect that scale-dependence decreases as more terms of the perturbative 
expansion are summed, as long as the series is converging. Is this realized by our results? In 
fixed-order perturbation theory, it is a common practice to vary the scale p, say, by factor two, 
and examine the stability of the prediction. It may be more natural to vary p such that £ changes 
by order y/T/N, as we argued at the end of Sec. 3.1. In either case, if we fix A^ in eq. (21), the 
variation in £ vanishes in the large A^ limit. A closer examination shows that the variation of 
t>c(p, N) — B(N), corresponding to these changes of p, also vanishes in the large A^ limit, as long 
as £ is close to 1. In this sense, our prediction becomes stable against scale variation at large 
orders. 

There exists an argument that the linear potential cannot emerge in perturbative QCD: From 
dimensional analysis, the coefficient of a linear potential should be non-analytic in as, i.e. of 



19 



order (A^ op ) 2 = p? exp[— An / (P as)]; therefore, it should vanish at any order of perturbative 
expansion. Within our large-order analysis, this argument is circumvented as follows. rV/v(r) 
includes terms of the form T n = { /3 ° a 2 s J^ log(/ir)} n for < n < N. If we substitute the relation 
(21) and take the limit JV -> oo while fixing n/N finite, it is easy to see that T n -> (A^ op r) 3 ^ N . 

Thus, perturbative terms converge to (A^ op r) p with positive powers < P < 3£. In fact, the 
power P has a continuous distribution. Our result shows that the continuous distribution can be 
decomposed into a sum of {r° , r l , r 2 , r 3 ^} terms, up to logarithmic corrections (for 2/3 < £ < 1). 
Non-analyticity in as enters through the relation (21). 

Thus, the characteristic feature of our large-order analysis is the prescription eq. (21). We 
may consider that an additional input has been incorporated through the relation (21) beyond 
a simple large-order analysis within perturbative QCD. Here, we emphasize that the number of 
parameters has not decreased from that of the original perturbative expansion (as, p, r and N) 
apart from N . (We fix A^jg, r and £ finite when sending the truncation order A" — > oo.) The 
"Coulomb" +linear part, t>c(p) + Cp, emerges independently of £ and A" in this limit. In this 
sense, we consider vc(p) + Cp a genuine prediction of perturbative QCD at large orders, within 
our estimate of the higher-order terms. 



3.3 Vn(v) at large orders: RG estimates 

In this subsection we examine V/v(r) for large A" using RG estimates of the all-order terms of 
Vqcd(^)- We examine three cases, corresponding to the estimates of Oy T (q) up to LL, NLL and 
NNLL in eqs. (5) and (6) [note that a n = P n (0) for n < 2]: 

(a) [LL] /3 , a : exact values, (3 n = P n (0) = (n > 1); 

(b) [NLL] (3q, Pi, a , a\\ exact values, (3 n = P n (0) = (n > 2); 

(c) [NNLL] p , Pi, Pi, ao, a±, ai. exact values, p n = P n (0) = {n > 3). 

Namely, cases (a),(b),(c), respectively, correspond to taking the sum up to n — 0,1,2 in eq. (5) 
and reexpanding in as{p)- From naive power counting of logarithms, one should also include US 
logarithms at NNLL. We examine them separately: 

(c') [NNLL'] Resummation of US logs is included via eq. (7), in addition to (c). 

We assume f3 , Pi, P2,a ,ai,a 2 (exact) > 0.* In the standard 1-, 2-, and 3-loop RG improvement 
of the QCD potential (in MS scheme), the same all-order terms as above are resummed; the 
difference of our treatment is that the perturbative series are truncated at 0(as(p) N )- We note 
that the estimate of higher-order behavior based on renormalon dominance hypothesis, as given 
in Sec. 2.3, is consistent with the above estimates, or more generally, with the RG analysis [38]. 
All the results for case (a) can be obtained from the results of the large-/3o approximation given 
in the previous subsection, if we replace A by A^°° p . 

Mb 

Below we summarize our results. (See Sec. 3.4 for derivation.) Similarly to the previous 
subsection, we can decompose V/v(r) into four parts: 

V N (r) = V c (r) + B(N, £)+Cr + V(r, N, f ) + (terms that vanish as N -> oo), (36) 
*This is the case when the number of active quark flavors is less than 6 and all the quarks are massless. 
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Figure 10: Integral contours C\ and C2 on the complex q-plane. q* denotes the Landau singularity of 
as{q). For 1-loop running, q* is a pole; for 2- and 3-loop running, q* is a branch point. In the latter case, 
branch cut is on the real axis starting from q* to —00. 



where 



V c (r) 



4irC F 2C F 



for 



Im 



7T 



dq a v (q), 

Ci W 



B(N,0 = lim ^Re / dqe* \a P v T (q) - [o^(q)] N }, 

r^O IX J Ci I J 

C P f prp , . 

C= 2^J C2 dqqav 

V(r, N, = ^v(r) - [V c (r) + + C r}. 



(37) 

(38) 

(39) 
(40) 



The integral contours C\ and C2 on the complex g-plane are displayed in Figs. 10(i),(ii). From the 
above equations, one can see that the "Coulomb" and linear parts, Vc(r) and Cr, are independent 
of £ and N, since Oy r (q) is independent of £ and N . 

The asymptotic behaviors of Vc(r) for r — > are same as those of Vqcd(^) in the respective 
cases, as determined by RG equations; the asymptotic behaviors of Vc(r) for r — > 00 are given by 
the first term of eq. (37) in all the cases. Namely, 



Voir) 
V c {r) 



2nC f 



1 



fo r|log(A M gr) 
4nC F 



1 - 



5 log|log(A M gr)| 



MA- 



ms 



for ' 



0, 



00, 



(41) 
(42) 



where 5 = in case (a). In the intermediate region both asymptotic forms are smoothly interpo- 
lated. 

Evaluating the integral eq. (39), the coefficient of the linear potential can be expressed ana- 
lytically in cases (a)-(c): 



C (a) = 
C (b) = 



2nC F 

~foT 
2nC F 

~foT 



V MS J 



A 2-loop\ 2 e 
MS J T(l + S) 



1 • ^6 ■«■■* 
fo 



e s 7(1 + 5, 5) 



(43) 
(44) 



where 7(2;, t) = JJ" (it t^ -1 e~* represents the incomplete gamma function; see e.g. [44] for defini- 
tions of A^? op . In case (c), the expression for C is lengthy and is given in App. C. Numerical 
values of C/A^ for various n ; are shown in Tab. I. 
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ni 


1 2 3 4 5 


cm/(a£">) 


2 


U./oz U.oii U.otw U.yol l.UUo l.Uyo 




2 


0.591 0.622 0.664 0.722 0.807 0.935 




■1 


1.261 1.317 1.385 1.465 1.556 1.644 



Table 1: Coefficients of the linear potential normalized by the Lambda parameter in MS scheme, 
for different values of ni. 

B(N,£) and V(r,N,£) depend on £ and diverge as N — > oo if £ is sufficiently large. In fact, 
apart from the overall normalization (and some details), behaviors of B(N,£) and T>(r,N,£) are 
similar to those presented in the previous subsection. We give two examples. 



V(r, JV, £) in case (b) with a x = and £ = 1: 



AC 



£=l,ai=0 /3o 



- rA?de opN i 3 r 2 

v MS ; 



1 / 3 \3<5/2 



1 / 3 VW2 



>J r(i + §5) 
f> °° r r(-x) 



12V2e/ r(l + §5) 
(log 2 + 7£ ) 



logiV + rf( b )(rAj± oop ) 



(45) 



3 \ 35 / 2 1 



2x(3-x) V2e7 r(l + ; 

(46) 



The asymptotic forms of S h \p) are given by 

1 / 3 \ 35/2 1 



[2 log | log p| +log| + 7 £ ] , p^0 or p -> oo, (47) 



12V2e7 r(l + §(J) 
and in the intermediate region both asymptotic forms are smoothly interpolated. 
• X>(r, oo,£) in case (b) with ai = and 2/3 < £ < 1: 



4C 



ai=0 



A) 



- CA^°P > i 3 S T .3f-i 

V MS ' 



x Re 



7o r(i + § 



5) 



r A 



2-loop\ * S ~mx/2 
MS y 



x5/2 



.(48) 



a;=3£— i s 



Its asymptotic forms are given by 

/^2-loop^ ^-1 



MS 



w 4C F r(-3Q . 
x — — — sin 



oi=0 



log(rAg op ) A> r(l + p) 



(?*) (I) 



3£<5/2 



\2eJ 

r — > or r — > oo. (49) 



The expressions when ai 7^ or in case (c) are more complicated and lengthy. 
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+ 



N = 100 
= 30 
N = 10 




A 2-loop 
MS 

Figure 11: [Case (b): NLL] V)v(r) for iV = 10, 30, 100 and £ = 1 (dashed lines). For comparison, 
the "Coulomb" +linear potential Vc(r) +Cr is also plotted (solid black). Constants have been added to 
V/v(r) and Vc(r) + Cr to make them coincide at r A^j!? OP = 0.5. We set n\ = 0. 



+ 




f = 0.9 



0.5 



1.5 



2.5 



r A 



2-loop 
'MS 



Figure 12: [Case (b): NLL] Vzv(r) for different values of £ and AT. (Dashed lines for £ = 0.9 and dot- 
dahed line for £ = 1.1.) For comparison, the "Coulomb" +linear potential Vc{r) +Cr is also shown (solid 
line). Other conventions are same as in Fig. 11. 
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0.5 1 1.5 2 2.5 3 

a 3-loop 
MS 

Figure 13: [Case (c): NNLL] V/v( r ) f° r N = 10) 30, 100 and £ = 1 (dashed lines). For comparison, the 
"Coulomb" +linear potential Vc{r) +Cr is also plotted (solid black). Other conventions are same as in 
Fig. 11. 




0.5 1 1.5 2 2.5 3 

r A 3 ' loop 
MS 



Figure 14: [Case (c): NNLL] V/v(r) for different values of £ and N. (Dashed lines for £ = 0.9 and 
dot-dahed line for £ = 1.1.) For comparison, the "Coulomb" +linear potential Vc{r) +Cr is also shown 
(solid line). Other conventions are same as in Fig. 11. 
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In Figs. 11-14, we show Vjv(r) for different values of A and £ in cases (b) and (c). (We set 
rii = in these figures.) They are compared with the "Coulomb" +linear potential V c (r)+Cr. The 
corresponding figures in case (a) can be obtained by simple rescaling of Figs. 6 and 8. Apart from 
the overall normalization, we see similar general features. Most importantly, Vjy(r) approximates 
well Vc(r) + C r at r < for a reasonably wide range of £ and A. For fixed £, Vjv(r) becomes 
steeper at r > as A increases [cf. eqs. (30), (45)]. For fixed A, Vn(t) is steeper for larger 
£ at r > this is because, if as(fj) is kept fixed and the truncation order is increased, all 

the higher-order terms additionally included contribute with positive sign. An only qualitative 
difference between case (a) and cases (b),(c) is that, for the same value of £ and A, Vjv(r) is slightly 
steeper (in comparison to Vc(r) +Cr) at r > in cases (b),(c) than in case (a). We postpone 
comparisons between cases (a),(b),(c) or comparisons with lattice computations of Vqcd(^) until 
Sec. 5. 

The effects of US logs in case (c') are very small (if we ignore shifts by r-independent constant). 
For instance, if we superimpose plots of Vjv(r) and Vc(r) + Cr of case (c') on Fig. 13, as we vary 
/ijA^ op between 2-5, they are hardly distinguishable from the corresponding lines of case (c). 

(Only Vjv(r) for A = 100 is visibly raised at rA^ op > 2.) The smallness of contributions from 
US logs stems from the small coefficient C\/(Q/3q) and suppression by \og[a s (q) /cxsifAf)] m eq. (7). 

Conclusions are essentially the same as those in the large-/3 approximation, because qualitative 
behaviors of V N (r) are similar: The "Coulomb" +linear potential, V c (r) +Cr, can be regarded as a 
genuine part of our prediction, while we may associate X>(r, A, £) with an 0(A^ CD r 2 ) uncertainty 
(and beyond) due to IR renormalons. Taking the variations of Vjv(r), corresponding to the different 
values of £ and A shown in Figs. 11-14, as a measure of uncertainties of the predictions for Vjv(r), 
the uncertainties are fairly small in the distance region r < A^. 

Let us compare our results in Sees. 3.2 and 3.3 with the results of the existing literature. 
The scale-fixing prescription according to the principle of minimal sensitivity was advocated and 
studied originally in [35]. In [45], a scale-fixing prescription close to eq. (21) was advocated, based 
on an analysis of large-order behavior of perturbative series a la renormalons; the prescription 
was used to suppress an ambiguity induced by UV renormalon, which is located closer to the 
origin than IR renormalons in the Borel plane. We studied in [30] the large-order behaviors of 
V N (r) using the scale-fixing condition eq. (21) but restricting to the case £ = 1, in the large- 
Po approximation and using the estimates by RG. Ref. [39] extended these analyses: Within 
the large-/3o approximation, and using the scale-fixing condition eq. (21), a general formula for 
the large-order behavior of a wide class of perturbative series was obtained and the relation 
to the Borel summation was elucidated; furthermore, the relation to the principle of minimal 
sensitivity was studied. Our present analysis is a direct extension of [30]; our results in the large- 
Po approximation in Sec. 3.2 are consistent with the general formula of [39] when the formula is 
applied to Vp (r). (Since the assumed singularity structure in the Borel plane is slightly different 
from that of Vp (r), slight modification of the formula is necessary). Unique aspects of [30] and 
the present work, besides being a dedicated examination of the QCD potential, are (a) the specific 
way of decomposition (close to Laurent expansion in r), and (b) inclusion of 2- loop and 3- loop 
running of as(q)- Furthermore, the separation into scale-independent (prescription-independent) 
part and scale-dependent (prescription-dependent) part is unique to the present work. 
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g/A 

Figure 15: (3q [a^ (q)]n an d Po&ysAq) y s- g/A for different values of £ and N. 



3.4 [ay T {q)]N as AT — ► oo 

In order to understand the properties of V/v(r) given in the previous two subsections, we examine 
behaviors of the truncated ^/-scheme coupling at iV — > oo, defined by 

[c£ T (g)]oo = lim [ay T (q)} N . (50) 

N^QO 

The relation (21) between «s(p) and N is understood in taking the limit. 
In the large-/?o approximation, one easily finds 



[ a vjh(q)]°o = Jim 



N^oo 



a s {p) 



1 — L 



N 



1 — L N 
lim a s (fj) r 



& r i-f4fi. (51) 



A)log(g/A) I ^ 

where L = ^|^log^^^j = 1 + ^logf^V There is no singularity at g = A, and [«v5 (g)]oo 
is regular at < |g| < oo in the complex g plane. For £ ~ 1, the first term in the curly bracket is 
dominant at UV (g — > oo), whereas the second term is dominant at IR (g — > 0). Hence, [a^g (g)]oo 
can be regarded as a modified coupling, regularized in the IR region, |g| < A; by including a 
power correction, the Landau pole of the original coupling, a{/^ (g) = 27r/[/3 log(g/A)], has been 
removed. (See Fig. 15.) 

One can test sensitivity of the prediction to the IR behavior of the regularized coupling by 
varying £. The results of Sec. 3.2 show that B(N) and D(p, N), associated with IR renormalons, 
are sensitive to the IR behavior of [«y^ (g)]oo, in accord with our expectation. On the other hand, 
^-independence of t>c(p) + Cp shows that t>c(p) + C p is determined only by the original coupling 
ocyfy(q), and that it is insensitive to the IR behavior of the regularized coupling. 

If we estimate the higher-order terms using RG, [a{/ T (g)]oo in case (a) is obtained simply by 
replacing A by A^ op in eq. (51). [«y T (g)]oo in case (b) can be analyzed using a one-parameter 
integral representation. To this end, let us first analyze the two-loop running coupling constant 
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as(q) (defined by eq. (6) with (3 n = for n > 2) and [a,s(g)]oo- Motivated by eq. (51), we separate 
[ q; s(q , )]oo into «s(g) and Aa 5 (g) = [a^g)]^ — as(g). They can be expressed in one-parameter 
integral forms, respectively, as 

/»oo /»oo 

®s(q) = / dxf!(x;q), Aa s (q) = - dxf^q), (52) 

JO J3£ 

with 

(See App. B for derivation.) These expressions are valid for a complex argument q if \q\ > g* = 
5 -5 / 2 A^ op ; they can be analytically continued to other regions by deforming the integral contour 
of x. as(q) and Aas(q), respectively, are singular at the Landau singularity g = g* (branch point), 
whereas their sum [as(g)]oo is regular at < \q\ < oo. One may find asymptotic behaviors of cts(q) 
and Aas{q) from the above expressions. The well-known asymptotic behavior of «s(g) as q — > oo 
is reproduced by rescaling xlog(g/A^ op ) — > x and expanding the integrand in 1/ log(g/A^°° p ). 
The asymptotic behavior of cts(q) as q — > 0, when g is varied along the path Ci of Fig. 10, can 
be obtained as follows. First rotate the integral contour clockwise around the origin, x = e~ ln y 
(0 < y < oo), then rescale y \og(A 2 ^ op / q) — > y and expand the integrand in 1/ log(A^ op /g). 
The asymptotic behavior of [a^g)],^ as q — > oo can be obtained by expanding the integrand of 
eq. (52) about x = 3£ except for the factor (g/ A^ op )~ x . The asymptotic behavior of [as(?)]oo as 
q — > 0, when g is varied along the path C±, can be obtained similarly, by first rotating the integral 
contour clockwise around the origin. The results read 

asiq) ~ Ml^-y (54) 

Aa s{q) (^ x gg, 9 ^ or 9 ^ oo. (55) 

SW A log( 9 /A|L~") \ q ) r(l + 3«i/2)' ; 

Thus, apart from the overall normalization, the leading asymptotic behaviors are identical with 
the one-loop running case, eq. (51). t 

ay T (g) in case (b) is given by ag{q) + (|^) a,s(g) 2 in terms of the two-loop running coupling 
constant. We also separate [a,s(g) 2 ]oo into a s (q) 2 and Aa s (q) 2 = [as(g) 2 ]oo — «s(<?) 2 , which can 
be analyzed similarly using the one-parameter integral expressions 

/»oo /»oo 

«s(<?) 2 = / dxf 2 (x;q), Aa s {q) 2 = - dxf 2 (x;q), (56) 

JO </3£ 



with 



f rr .^ _ 87r2 ^/A^ioopA^ r ^/2 7(1 + ^/2,3^/2) 

^ r(i + x5/2) • (57) 



The asymptotic forms of Aats(q) are given by 



2^ ( tfjjT\ 3? v 4VTA, ^ H5/ 2 7(1 + 3^/2,3^/2) 

' ' HS/2) ' 

or g — > oo, (58) 



Aa * (g)2 ~ ftMC P ) (, 7 J x A " r(i + 3^/2) 



t There are qualitative differences in the subleading asymptotic behaviors. 



27 



while the asymptotic behaviors of as(q) 2 are given simply by the square of eq. (54). Thus, 
(f 1 ) a s(q) 2 term does not affect the asymptotic behaviors of ay T (q) and Aa P 7(q) = [ay T (<?)]oo — 
Oi v T {q)-, apart from the overall normalization. 

The same is true in case (c). The leading asymptotic behaviors of a v T (q) and Aa v T (q) as 
q — > (when q is varied along the path C\) and q — > oo are same as those in case (a) (one-loop 
running), besides the overall normalization. These IR and UV behaviors determine the behaviors 
of Vjv(r) at r — > oo and r — > for large N. For this reason, the truncated potentials Vjv(r) 
have qualitatively similar features in all the cases examined in Sees. 3.2 and 3.3. We also note 
that, since [tty T (5')]oo has no singularity along the positive real axis, and since Aa v T (q) is more 
dominant than a v T (q) in IR, the leading IR behavior of [tty T (5')]oo, when q is sent to +0 along 
the positive real axis, is same as that of Aa v T (q) as q — > (when q is varied along the path Ci). 



3.5 How to decompose V/v(r) 

We explain how we decompose Vjv(r) into 4 parts, as given in Sees. 3.2 and 3.3. Integrating over 
the angular variables in eq. (20), one obtains 

2C F f°° , singr r PT , 2C F T f°° , e iqr . PT/ 

V N {r) = / dq -[a v T (q )} N = Im / dq [a v T (q)] N . 59 

7T Jo qr vr Jo qr 

We separate the integral into two parts, according to the different asymptotic behaviors of 
["y T (<?)]oc, i.e. a v T (q) and Aa v T (q) = [a{7(<?)]oc - a v T (q): 

V N (r) = U 1 (r) + U 2 (r,N,Z), (60) 

U 1 (r) = -^lm / dq — a v T (q), (61) 
J Cl QT 

U 2 (r,N,0 = -^lm j c dq^—{[a v T (q)] N -a v T (q)}. (62) 

We deformed the integral contour in order to avoid the Landau singularity on the positive real 
axis; see Fig. 10(i). Contributions from the Landau singularity cancel between Ui and U 2 , since 
the original integral (59) does not contain the singularity. 

Since [ay T (q)] N — ay T (q) ~ ? -3 V logg as iV — » oo, the integral in eq. (62) becomes IR divergent 
in this limit for £ ~ 1. On the other hand, the negative power of q induces the positive power 
behavior of r in U2 in the large N limit. Assuming £ > 2/3, let us define 

A 

U 2 (r, N,g) = — + B(N, + Cr + V(r, N, f) + (terms that vanish as N -> 00), (63) 
r 

where V(r,N,£) is subleading as compared to Cr at short-distances. A and C can be extracted 
as follows. 

A — lim r U 2 — lim Im / dq \ [a v T (q)] N — a v T (q) \ 

r^O r-*0 71 q I J 

= / ^{[^*-»?( 5 )} = -^ / *^ =-^, (64) 

™ Jc 2 V 1 > 1X1 JC 2 1 A) 

C = ^ n ~(rU 2 ) = ^m^lm f dqe^q \[a v T (q)] N - a v T (q)\ 

r-»0 Z Or r^O 7T Jq I J 

= I dqqa v T (q). (65) 
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To show the last equality of eq. (64), we may use the RG equation (6), or, we can evaluate the 
integral explicitly using a v T (q) at LL, NLL and NNLL. [a:y T (g)]Ar does not contribute because 
it has no singularity inside the contour C 2 , hence, both A and C are independent of £ and N. 
Similarly, B and V are given by 

£(iV,0 = lim|-(rf/ 2 ) = lim-^Im ? / dqe* \[a v T (q)} N - a? T (q)\, (66) 



V(r,N,0 = U 2 (r,N,0 

2C F T f , e iqr - [1 + iqr + \(iqr) 



- + B(N,0+Cr 
r 



J Cl QT I J 

In the limit iV — > oo, B is IR divergent, and if £ > 1, X> is also IR divergent. We would want to 
factor out divergent part as iV — > oo in these cases.''" 

If £ < 1, T> is finite as iV — > oo. Then, we may insert the expression for Aay T (q) obtained in 
the previous subsection. In case (a) or in the large-/3o approximation, it is convenient to deform 
the integral contour into the upper half plane by setting q = ix/r (0 < x < oo). Then, one readily 
obtains eq. (33). To find the asymptotic forms, eq. (34) and subleading terms, one expands the 
integrand (inside Im[...]) by log a;. In cases (b) and (c), we may insert the integral expressions for 
Aay T (q) to eq. (67) and integrate over q. Thus, we obtain one-parameter integral expressions for 
V, except for the coefficient of a 2 in case (c).§ We may find asymptotic forms of V, for instance, if 
we insert the asymptotic expansion of Aay T (q) [eqs. (55), (58) and subleading terms] and proceed 
as in case (a). 

When £ = 1, we may factor out the divergence as 

7T J a Q r 5=1 

-^Im/^ g M( [a PT (?) ],_ a PT (g) } (68) 
7r J 6qr I ) £=i 

where go is an IR cutoff to remove the IR divergence as q — > 0. In all the cases (a)-(c), one may 
extract the divergent part from the second term, which may be taken as 

f 1 a 2 r / 3 \ N i 1 f 1 \ 

L^K 1 -*** 9 ) -i]— (l0g/v + l0g2+%) + o (7F)' (69) 

apart from an overall normalization (proportional to r 2 ). Here, we have rescaled q to a dimen- 
sionless variable q. If £ > 1, one may factor out the divergences in a similar manner; one should 
subtract powers of r as many times as needed to remove all the IR divergences. It is even simpler 
to factor out divergences from B(N,£); for instance, see eq. (21) of [30] for £ = 1 and in the 
large-/3o approximation. 

We define the "Coulomb" potential (with logarithmic corrections at short-distances) as 

V c (r) = U 1 (r) + j. (70) 

It is determined by ay T (q), therefore, it is independent of £ and N. Since the leading behavior 
of Vjy(r) as r — > is const. /(r log r) as determined by RG equation, the A/r term of U 2 must 

*The integrals (64)-(67) are convergent in the UV region, assuming the double limits eq. (23). 
§We were able to reduce the coefficient of a,2 in case (c) only to a 2-dimensional integral form. 
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be cancelled by the 1/r term contained in U\. To compute the asymptotic forms of U\(r) as 
r — > and r — > oo, one may insert the integral expressions for ay T (q), given in the previous 
subsection, and integrate over q; then rescale x log(l/rA^ op ) — * x and expand the integrand in 

1/ log(l/rA^ op ). In this method, however, one should carefully choose the integral contour for x 
to avoid singularities. Another method is as follows. Consider first the case (a). We deform the 
integral contour into the upper half plane on the complex g-plane and integrate by parts: 



AC F , f°° dx e~ x 

x log(x/ p) + in/2 



AC 



F 



Im J dxe x log (^log p — log x — —J, (71) 



where p = rAi^ op . By expanding the integrand in logx, we obtain the asymptotic forms eq. (27). 
In cases (b) and (c), we may proceed in parallel with the above steps, after appropriate change of 
variables in the integration. Then eqs. (41) and (42) are obtained. 

Let us summarize our algorithm for decomposing Vn(t) into Vc(r) + B + Cr + V(r). First we 
separate [«{/ T (5 , )]at into 2 parts according to the different asymptotic behaviors of [«y T (g , )]oo as 
q — > and q — > oo. The separation is particularly simple in the one- loop running case, eq. (51). 
Next we deform the integral contour into upper-half plane to avoid the Landau singularity. Thus, 
V/v(r) is separated into U\{r) and [^(r). ^(r) has a power series expansion in r from 0{r~ v ) 
to 0(r). Beyond that order, power series expansion in r breaks down due to non-analyticity in 
r.^ Hence, U 2 (r) is naturally decomposed into A/r + B + Cr + V(r). Then, A/r is combined 
with Ui(r) to form Vc(r), which behaves as a Coulombic potential in the entire range of r apart 
from logarithmic corrections. Thus, V^{r) is decomposed into {Vc{r), B,Cr,T>(r)}, which corre- 
spond to {r _1 , r°, r 1 , r 3 ^ -1 } terms; the r~ l and r 3 ^ -1 terms include logarithmic corrections. The 
"Coulomb" +linear potential, Vc(r) +Cr, is determined only by «y T (g), hence, it is independent 
of f and N. 

If we choose a different prescription to avoid the Landau singularity in defining U\{r) and 
U2(r), each of them will change (but the sum V)v(r) will not). Consequently, Vc(r) and T>(r) will 
also change. Since, however, the contribution of Landau singularity does not have a simple power- 
like form in r, in other prescriptions Vc(r) = Ui(r) + A/r is not Coulomb-like at large-distances 
(rather oscillatory). We consider our decomposition natural, in the sense that it is closest to a 
decomposition into terms with simple powers in r (cf. argument at the end of Sec. 2.4), as well as 
since Vc(r) + Cr is a good approximation of VJv(r) at r < A^L for a reasonably wide range of £ 
and N (cf. Figs. 8,11-14). 



4 Renormalization Schemes in OPE of Vqcd(^) 

We analyze Vqcd(?") in OPE when r <^ Aq^ d . As we have seen in Sec. 2.4, the leading short- 
distance contribution is given by the singlet potential Vg(r). In this section, we provide renor- 
malization prescriptions for Vg(r) explicitly. We also show that the 'Coulomb" +linear potential, 
extracted in the previous section, can be qualified as a singlet potential (Wilson coefficient) de- 
fined in a specific renormalization scheme. These renormalized singlet potentials are free from 
IR renormalons and IR divergences; hence, they can be computed systematically and (in princi- 
ple) we can improve the predictions to arbitrary precision. Correspondingly, the non-perturbative 
contributions are unambiguously defined. 

^ For simplicity, we restrict ourselves to the case 2/3<£<lin this and the next paragraph. 
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4.1 Factorization scheme vs. "Coulomb" + linear potential 



Following the argument of Sec. 2.4, let us define a renormalized singlet potential, in a scheme 
where the IR divergences and IR renormalons are subtracted, as 



2C f 



7T 



dq 



sin(gr) (R) 



a 



with 



a 



(R) 
V S 



a 



qr 



v s 



Say{q;fXf). 



(72) 



(73) 



5ay(q; fJ>f) is the counter term which subtracts the IR divergences of a v (q), given as multiple 



poles in e. (We assume that a v (q) is computed in dimensional regularization.) 

Let us consider two schemes in particular for defining Say(q; /if). First one is to subtract the 
IR divergences of a v T (q) in the MS scheme. Explicitly, at NNNLO, we set 



5a v (q;n f ) = a s (^) 



47T 



~ - 8 log Q + A lE - 4 log(47r) + 2 log 



, (74) 



where we retained (only) the physical US logarithm according to the argument given below 
eq. (119). Here, /i/ represents the scale at which loop momenta are effectively cut off. The loga- 
rithms induced by running of a$ can be resummed up to NNNLL by setting // — > q in ay^(q; jj,f). 
As we saw in Sec. 3.3, resummation of US logarithms does not give sizable effects, so we will not 
try to resum US logs but rather include US logs only up to NNNLO, as given in eq. (74).* 

Second scheme is to regularize the IR divergences by expanding a v T (q) as a double series in 
as and log as- Then, no artificial subtraction from the IR region of loop momenta is made, 
and aty9(q;iif) becomes independent of /i/. At NNNLO, 5ay(q;j2f) is obtained from the Fourier 
transform of eq. (120): 



6a v (q;n f ) = a s {n) 



x 72vr 2 



X - - 4 |21og(^) + log(47r)| + 6 7s - \ + 2 log^o^/i)) 



• (75) 



Indeed it is independent of ///. This prescription introduces a physical scale AV(r) = Vo(r) — 
Vs(r) ~ CaolsI t as an IR regulator in loop integrals, hence, contributions from q < AV(r) are 
suppressed [1]. Below, we will resum powers of log(/i/g) associated with the running of as but 
not powers of US log 0:5, for the same reason as in the first scheme. In Sec. 5.1, we will compare 
the two schemes eqs. (74) and (75) numerically. 

Dependence of vJ R \r; fij) on fif is introduced through subtraction of the IR divergences (in the 
first scheme) and of the IR renormalons. The subtraction of the IR divergences in the first scheme 
induces logarithmic dependences on /j,f [eq. (74)], while the subtraction of the IR renormalons 
induces power-like dependences on /if. The former resides in the counter term 5av(q; (J>f), and 
the latter arises from the lower cutoff of the integral in eq. (72). Roughly, 



d 



(76) 



* Resummation of US logs up to NNLL is achieved if we omit log(/it//g) in eq. (74) and make the replacement 
eq. (7). Resummation of US logs up to NNNLL has not been computed yet. 
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neglecting the r-independent part. Note that /i/r <C 1 due to the hierarchy (15). 

Corresponding to the above definitions of Vg R \r; ///), we define the "Coulomb" +linear potential 
from a { yj(q;/j, f ) by 

V c+h (r) = V^\r)+C^r, (77) 

where 

Wr)~& / c/9 ^-^. m / Ci ^W,), (78, 

C (R) = ^^g«gW/)- (79) 

Up to NNLL, it is natural to take 5ay(q] A 4 /) = 0> hence, in this case, Vc+h( r ) coincides with the 
"Coulomb" +linear potential obtained in the large-order analysis, eqs. (37) and (39) [cf. eq. (64)]; 
in particular, Vc+h( r ) is independent of /if up to this order. 
Below we will show that 

KfV;^/) - V c+ L(r) = const. + 0(fi 3 f r 2 ). (80) 

Eqs. (17) and (80) imply that, within the framework of OPE, short- distance contributions (q > /!/) 
determine the "Coulomb" +linear part of the QCD potential, hence it is predictable in perturbative 
QCD. The residual term (apart from an r-independent constant) is of order /i^r 2 , which mixes 
with the non-perturbative contribution SE\js(r), and is subleading at r <C /ij 1 . These features 
are consistent with our expectation discussed at the end of Sec. 2.4. 
Eq. (80) can be shown as follows. According to eqs. (72), (77)-(79), 

VfVi W )-V c+L (r) 

= c, r ^W) + 2c, Im i ag) _ C(R) 

mr Jc 2 q ^ Jc 3 1 r 

where the integral path C3 is shown in Fig. 16. Since Hfr <C 1, we may expand the Fourier factor 
as e iqr = 1 + iqr — \{qr) 2 + ... in the integral along C3. Then the leading term of the expansion 
cancels against the first term of eq. (81), while the third term of the expansion [— \{qr) 2 } cancels 
against — r of eq. (81). Therefore, only remaining terms on the right-hand-side of eq. (81) are 
const. + C>(/r}r 2 ). 

One may think that in defining Vg (r;///) subtracting the integral eq. (18) is not sufficient 
for subtracting all the IR renormalons. The relation eq. (80) is unchanged, even if one subtracts 
the IR renormalon contributions using whatever other sophisticated method for estimating them. 
This is because the IR renormalons in Vg R \r; /ij) take the form const. + (9(AQ CD r 2 ). 

The perturbative expansion of Vg K \r;/j,f) may still be an asymptotic series (due to e.g. UV 
renormalons^). Since the IR renormalons have been subtracted and the factorization scale is set 
as fif 3> Aqcd, we may expect that Vg R \r; fif) is Borel summable.* (At least, the Borel integral is 
convergent in the large-/9 approximation.) Then, we may define Vg R \r; /ij) from the perturbative 
series either by Borel summation or according to the prescription of [45, 39]. Thus, Vg R \r; ///) 
can be computed systematically (based on perturbative QCD). 

t Nevertheless, we note that up to now UV renormalons have not been identified in Vqcd(^)- 
1-This is up to the uncertainties caused by the instanton-induced singularities in the Borel plane, which we 
neglect in our analysis. 
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Figure 16: Integral path C3 in the complex g-plane. g* denotes the Landau singularity of as(o)- 
For 1-loop running, is a pole; beyond 1-loop running, g* is a branch point. In the latter case, 
branch cut is on the real axis starting from g* to —00. 



4.2 Vc+l(^) as a //^—independent renormalized singlet potential 

Up to NNLL, the "Coulomb" +linear potential Vc+l(?") = Vc( r ) + Cr was extracted from Vjv(r) as 
a prescription-independent part, corresponding to a renormalon-free part, in Sec. 3. We have also 
seen that V c+ ^(r) coincides, up to 0(r 2 ), with Vg R \r; ///), which is the Wilson coefficient free 
from IR renormalons. Therefore, unlike original perturbative expansion of Vqcd(?")j we expect 
that Vc+h(r) is free from intrinsic uncertainties. In fact, Vc+l(?") can be computed systematically 
and its accuracy can be improved (in principle) to arbitrary precision as follows. Since a^(q; (if) 

does not contain IR renormalons,* we may compute a^(q; /J,f) and improve accuracy of its value 
along the contours C\ and C 2 by including corrections at LL, NLL, NNLL, and so on. (We further 
improve the series by Borel summation or by the prescription of [45, 39] if necessary.) In the 
infinite-order limit, a^\q;fif) is expected to be finite everywhere along these contours."'' Then, 
via eqs. (78) and (79), Vc+l(?") can be computed. (It is beyond our scope to prove convergence of 
Vc+h(f)- Here, we have eliminated all of the known sources of divergences.) 

In the definition of Vc+h(f), renormalization scheme dependence enters through the definition 
of aty£(q;Hf) or of 5a>v(q; At/). Below, we will focus on the second scheme, discussed in the pre- 
vious subsection: To regularize the IR divergences by expanding ay T (q) as a double series in as 
and logas. Then, V c +l(^) becomes independent of /if. In view of the construction of OPE of the 
QCD potential, the "Coulomb" +linear potential Vc+h( r ) defined in this scheme can be qualified, 
without any problem, as a renormalized singlet potential (Wilson coefficient) defined in a specific 
scheme. We remind the reader that the bare singlet potential Vs(r) coincides with the perturba- 
tive expansion of Vqcd(^), and that the renormalized singlet potential is defined by subtracting 
IR renormalons and IR divergences from it. According to the large-order analysis and the def- 
inition (77)-(79), Vc+L( r ) matches this requirement. Furthermore, eq. (80) ensures consistency 
of identifying Vc+l(?") as a renormalized singlet potential. Vc+h{ r ) is well-defined, systematically 
computable, and free from ambiguities induced by IR renormalons or IR divergences. The only 

*In a more sophisticated estimate of renormalons than eq. (18), one may find that a Vs (q;fif) still contains 
IR renormalons. If this turns out to be the case, according to our philosophy, it is appropriate to subtract the 
renormalons by modifying the counter term 5av(q; A*/)- 

1 There is no known source of divergence except the instanton-induced singularities, which we neglect in this 
paper. 
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notable difference from the ordinary OPE is that, in this scheme, Vc+h(f) is independent of the 
factorization scale /if. 

Eq. (80) shows that Vg R \r; fj,f) approaches Vc+l(?") as we reduce /Uj(^> Aq CD ). This matches 
our naive expectation: Vc+i,(r) is obtained from the bare Vs(r) by subtracting the part corre- 
sponding to IR renormalons, which reside in the region q ~ Aqcd; on the other hand, Vg (r; /if) 
is obtained by cutting off a larger domain q</if. 

4.3 SEus(r): C+L scheme and factorization scheme 

We may replace Vg(r) in eq. (17) by Vc+h( r ) and Vg R \r; /if), respectively, and define the non- 
perturbative contributions 5E\js(r) corresponding to both schemes. Alternatively,-'' via eq. (16), 
we may identify them, respectively, with Vqcd(^) — Vc+l(?") and Vqcd(^) — Vg (r;///). Let us 
call the former as 5Eys( r ) in the C+L scheme and the latter as 5Eus(r) in the factorization 
scheme. As mentioned above, there are no intrinsic uncertainties in Vc+l(?") and V s (r;///), so 
that 5E\js(r) in both schemes are unambiguously defined. IR renormalons have been subtracted 
from the bare Vg(r) and absorbed into SEjjs(r). 5E\js(r) in the C+L scheme is independent of 
fif, while SEjjs(r) in the factorization scheme depends on fif. 

According to the argument in Sec. 2.4, r-dependence of 5E\js(r) in the factorization scheme 
can be predicted. (We always set /if 3> Aq C d in the factorization scheme.) It is either of order 
fijr 3 (very small r) or of order fi^r 2 (small r), depending on the relation between AV(r) and /if. 
We expect the latter r-dependence in the distance region of our interest. On the other hand, the 
r-dependence of 5E\js(r) in the C+L scheme can be estimated in parallel with that of 5E\js(r) in 
the factorization scheme with /if ~ Aqcd- Therefore, r-dependence in the C+L scheme can be 
predicted for small r (corresponding to AV(r) 3> Aqcd) to be order AQ CD r 3 . If r is not sufficiently 
small (AV(r) ~ Aq C d), precise r-dependence is not known. Since, however, SEu S (r) in the C+L 
scheme contains no other scale than Aqcd> it should be at most order Aqcd at r < Aq CD . Namely, 
it is much smaller than 5Ejjs(r) in the factorization scheme (order /x^r 2 ), provided ///is sufficiently 

large. It means that Vc+l(?") is much closer to Vqcd(^) than Vg R \r;///) in the distance region of 
our interest. This gives us a good motivation to analyze Vc+l(?") in OPE, in addition to the more 
conventional factorization (/ij-dependent) scheme. 

All the above arguments are based on order-of-magnitude estimates. We would like to make 
the statements clearer by making a quantitative analysis. 

5 Determinations of di^usM and r$ Aj^g 

In this section, we compare the singlet potentials, which we defined in different schemes in the 
previous section, and recent lattice data for the static QCD potential. According to previous 
analyses, (a) accuracy of the prediction for Vg R \r; [if) or Vc+h(f) can be improved systemati- 
cally; (b) the difference 5E vs (r) = Vq C d(0 - V^ R) (r;/i/) [Vq C d(0 - Vc +L (r)] is expected to be 
(9(/ijr 2 ) [(9(AQ CD r 3 ) at small distances, whereas precise form is unknown at larger distances] and 
is non-perturbative. We verify these properties numerically by comparison to lattice computa- 
tions of Vqcd(?")- Then we determine the size of the non-perturbative contribution 5E\js(r). As a 

*As for SE\js{r) in the factorization scheme, contributions from gluons close to the UV cutoff q ~ [if can be 
computed reliably in expansion in as using eq. (17) (although the entire SE\js(r) cannot be computed reliably). 
Then, one can show explicitly that the /iy-dependence of Vg R \r; fj,f) [cf. eq. (76)] are cancelled by the yUj-dependence 
of SEus(r) [24], showing consistency of defining <5_E us (r) in two ways. 
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byproduct, we determine the relation between A^g- and lattice scale (Sommer scale) at the same 
time. 

Throughout this section, we use lattice data in the quenched approximation, since in this case 
lattice data are most accurate in the short- distance region. In computations of Vg R \r; nf) and 
Vc+l(?") : we set ni = accordingly; except where stated otherwise, we take the input parameter 
as a s (Q) = 0.2, which corresponds to* A^°° p /Q = 0.057, Aj±°° p /Q = 0.13, Ajg op /Q = 0.12; 
at NNNLL, except where stated otherwise, we use the estimate of a 3 by Pineda in eq. (119), 
a 3 = 292 x 4 3 = 18688 [26]. An arbitrary r-independent constant has been added to each 
potential and each lattice data set to facilitate comparisons in the figures. Methods for numerically 
evaluating Vg R \r; /iy) and Vc+l(^) are shown in App. D. 

We relate the scale for each lattice data set to Aj^g in the following manner. For each lattice 
data set we calculate (or use the given value of) the Sommer scale r defined by 

r 



dr 



= 1.65. (82) 



(For reference to the real world, it is customary to interpret ro = 0.5 fm £s 2.5 GeV -1 .) Then the 

lattice data are expressed in units of ro. In Sec 5.1, we convert the units into A^ op using the 
central value of the relation 

r Ag op = 0.602 ± 0.048, (83) 

as obtained by [47]. In contrast, in Sec. 5.2, we will not use the relation between r and A|^ op 
as an input but rather determine this relation from a fit to the data for 5E\js(r). We will explain 
the mechanism why this is possible. 



5.1 Consistency checks 

Here, we verify various properties of the singlet potentials Vg R \r;/jLf) and Vc+l(?") an d of the 
corresponding non-perturbative contributions. 

First we compare the "Coulomb" +linear potential Vc+l(?") up to different orders, in the fif- 
independent scheme defined in Sees. 4.1 and 4.2. Up to NNLL, they coincide with Vc{r) +Cr of 
Sec. 3.3. We also compare them with lattice calculations of the QCD potential. See Fig. 17. * We 
see that Vc+h( r ) up to different orders agree well with one another at small distances, whereas at 
large distances V c+ l(^) becomes steeper as we include higher-order terms via RG; cf. Tab. 2. This 
feature is in accordance with the qualitative understanding within perturbative QCD, in which 
the potential becomes steeper due to the running of the strong coupling constant, since <yy T {q) 
increases more rapidly at IR as we include higher-order terms. The lattice data and Vc+h(f) 
also agree well at small distances, while they deviate at larger distances. More terms we include 
in V c+ L(r), up to larger distances the potential agrees with the lattice data.* Theoretically, we 
expect V c+ L(r) to converge as we increase the order. The current status seems to be consistent 

*As well known, when the strong coupling constant at some large scale, e.g. as(m&), is fixed, the values of 
^l-loop ^2-loop ^3-loop jjjggj. substantially. As a result, if we take a common value of A-ria as the input 

MS ' MS ' MS J ' Mb ^ 

parameter, Vc+lM up to different orders differ significantly at small distances, where the predictions are supposed 
to be more accurate. 

*If we use Chishtic-Elias's estimate of S3, the NNNLL line in Fig. 17 hardly changes. If we use the estimate of 
(Z3 by large-/3o approximation, the NNNLL line is located between the present NNNLL line and NNLL line. 

* It is worth noting that the NNLL line in Fig. 17 is numerically very close to the NNLO prediction obtained 
with principle of minimal sensitivity in [27]. 



35 



O i 

o 



2 


-2 
-4 
-6 



NNNLL (est.) 




0.2 



0.4 



0.6 

. A 3-loop 

' MS 



0.8 



Figure 17: Comparison of Vc+h( r ) in the ///-independent scheme (solid lines) and the lattice data in 
quenched approximation [Takahashi et al. [48] (o), Necco/Sommer [49] (•), and JLQCD [50] (*)]. Input 
parameters for Vc+h(f) ar e cxs(Q) = 0-2 and n\ = 0; at NNNLL, Pineda's estimate for 03 is used. 





LL 


NLL 


NNLL 


NNNLL 




0.1836 


0.6950 


1.261 


1.758 



Table 2: Coefficient of the linear potential [eq. (79)] in units of (A|^ op ) 2 . Conventions are same as in 
Fig. 17. 
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Figure 18: Comparison of lattice data and Vg (r; /U/) up to NNLL for ^j/A^ op = 2 to 5 (downwards). 

As a reference, we also plot Vc+l(t) up to NNLL. The lattice data and parameters for V s (r;^j) and 
Vc+l(^) are same as in Fig. 17. 

with this expectation, since the lines in Fig. 17 apparently converge to the lattice data. One may 
adjust the input value ag(Q) such that convergence becomes fast at some particular r; present 
choice as{Q) = 0.2 leads to a fast convergence at rA^ op < 0.1. We may increase the value of 
input ots{Q) such that convergence becomes fast at larger r. Then, Vc+h( r ) up to different orders 
come closer to one another at r A|^? op > 0.1. [The relation between Vc+l(?") up to NNLL and the 
lattice data remains unchanged, since we use the 3-loop RG relation to fix the lattice scale, i.e. 
this relation is invariant under 3-loop RG evolution.] 

Next we compare the renormalized singlet potential Vg (r; ///) for different values of the fac- 
torization scale /j,f. In Fig. 18, we plot Vg (r; /if) up to NNLL for /ij/A^ op = 2,3,4,5. Note that, 
up to NNLL, there is no distinction between the first scheme and the second scheme of Sec. 4.1. 
We see that Vc+l(V) is located closer to the lattice data than Vg (r; (if) for all ///. As we vary fif, 
the variation of Vg (r; /i/) is larger at larger distances. As we lower /!/, Vg (r;/i/) approaches 
Vc+h{f)- These features are in agreement with the argument given in Sec. 4.3. Note that we 
cannot lower jjf below the Landau singularity g^" loop ~ 1.53 A^ op . Vg (r;^t/) up to different 
orders behave similarly: As /// is lowered, Vg R \r;/j,f) is raised at larger distances, approaching 
Vc+h(r) up to the corresponding order. For fixed fif, Vg (r;///) agrees with lattice data up to 
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Figure 19: <5£ us (r)/Ag op vs. r Ag op . The lattice data from Tab. 2 of [49] are used. Parameters for 

Vc+l( t ) and (r;/j,f) are same as in Fig. 17. Lines represent fits to the data points in the range 
r A|^? op < 0.5 by third-order polynomials, given explicitly in Tab. 3. (a) C+L scheme, up to LL, NLL, 
NNLL and NNNLL; C+L scheme in the large-/?o approximation is also shown, (b) first scheme within the 
factorization scheme (subtraction of IR divergences by MS scheme); for display purposes, 5E\js(r)/ A?-^ op 
at NNLL are shifted by +1 vertically. 



larger distances as we increase the order. 

We turn to the measurement of 5Eu S (r). In computing 5Eu S (r), we use the lattice data from 
Table 2 of [49] for a non-perturbative computation of Vqcd(?")> since this data set seems to be most 
accurate at short distances. In Fig. 19(a), we plot 5E\js(r) in the C+L scheme [Vqcd(^) — Vc+l(?")] 
in units of A|^ op . The errors of the data points, due to the errors of the lattice data for Vqcd(?")> 
are comparable to or smaller than the sizes of the symbols used for the plot. Also shown in the 
same figure are fits to the data points of the form A\ p + Ai p 2 + A% p 3 , where p = rA^ op . Only 

the data points in the range r A|^ op < 0.5 were used for the fits. We have added r-independent 
constants such that all the fits go through the origin. As we increase the order, <5£\js( r ) becomes 
smaller. We see that the cubic fit becomes a better approximation in a wider range in r A^ op < 1 
as we increase the order. [Here, we determine 5E-[js(r) m expansion in r, hence, we fit the data 
in the small r region (r A|^ op < 0.5); it helps to enhance sensitivity to the coefficients of r n for 
small n. Since 5Eu S (r) ~ 0(AQ CD r 3 ) at sufficiently small r, i.e. 5Eu S (r) — > as r — » 0, it would 
make sense to perform a polynomial fit; however, in the next subsection, we reconsider this naive 
picture and give a more complete analysis.] 

We can learn more detailed features from the explicit polynomials obtained from the fits, 
shown in Tab. 3. The coefficient of the linear potential decreases as we increase the orders, from 
LL to NNNLL. Up to NNNLL, the coefficient in units of (A^° p ) 2 is about 0.7. We may compare 
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C+L scheme factorization scheme factorization scheme 

= 3A ^° P ) = 5A ^° P ) 


LL 

NLL 

NNLL 

NNNLL, 1st scheme 
NNNLL, 2nd scheme 
Large-/3o appr. 


6.7p-9.6p 2 + 8.7p 3 
3.5p-3.6p 2 + 3.6p 3 

1.6p-0.3p 2 + 0.8p 3 1.6p + 0.9p 2 + 0.8p 3 1.5 p + 3.5 p 2 - 0.8 p 3 

0.7p+ 1.8p 2 + 0.1p 3 0.6p + 5.0p 2 - 1.8p 3 
0.7p + 0.8p 2 -0.1p 3 0.7p + 2.0p 2 -0.0p 3 0.6 p + 5.1 p 2 - 1.9 p 3 
-2.7p+10.9p 2 -9.1p 3 



Table 3: Fits of 6E\js(r)/ A^ op by cubic polynomials in the region p < 0.5, where p = rA^ op . (See 

Fig. 19 for plots.) The lattice data [49] are used. Parameters for Vc+l(0 and Vg (r; p/) are same as in 
Fig. 17. 



this value with the string tension (coefficient of linear potential) extracted from the large-distance 
behavior of the lattice data a « 3.8 (A^- op ) 2 [49]. Thus, the linear potential in 5Eus(r) is quite 
small comparatively at our current best knowledge. We are interested in SE\js(r) in the infinite- 
order limit. Up to our current best knowledge, however, SE\js(r) in the C+L scheme has not 
stabilized yet, although we see a tendency that it approaches a quadratic form. We conclude 
that the present status is consistent with 5Eu S (r) in the infinite-order limit being order AQ CD r 2 at 

r A|^ op < 1 (vanishing linear potential). From Figs. 17 and 19(a), it is not clear whether the limit 
would also be consistent with order Aq CD r 3 or with zero. We will clarify this point quantitatively 
in the next subsection. 

In Fig. 19(b) are plotted 5Eu S (r) in the first scheme (IR divergences of a v T (q) are subtracted in 
MS scheme) within the factorization scheme [Vqcd(?") — Vg R ^(r; p^)]. In this figure and in Tab. 3, we 
see that SE-[j S (r; pj) approximate the quadratic form 0(p 3 r 2 ) expected from OPE. Approximation 
to the quadratic form is more evident than in the C+L scheme, since the coefficents of the quadratic 
terms are larger. Differences between the first and second scheme in the factorization scheme are 
tiny: If we plot, in the same figure, SE\js(r) in the second scheme (IR divergences of a v T (q) are 
regularized by double expansion in as and log as), they are hardly distinguishable from the lines 
of first scheme. This is also confirmed by the explicit forms of the polynomial fits in Tab. 3. 

We may compare the fits of 5Eus(r) in the factorization scheme and those in the C+L scheme 
in Tab. 3. The coefficients of linear terms are almost same in both schemes, up to NNLL and up to 
NNNLL, respectively. This confirms consistency with eq. (80). The difference of the coefficients of 
quadratic terms between both schemes is expected to be proportional to p 3 in the limit p/ ^> Aq CD , 
according to eq. (80). This relation is roughly satisfied in Tab. 3 as well. § 

Next we compare Vc+l(?") in the large-/3o approximation,^ analyzed in Sec. 3.2, with the lattice 
data. We take two different values of the input parameter: as(Q) = 0.2 and as(Q') — 0.5, which 
correspond to Q/A^ op = 17A and Q'/A^ op = 3.1, respectively. See Fig. 20. Since the large- 
f3o approximation incorporates only 1-loop running of as, the prediction for Vc+h( r ) is fairly 
scale dependent." Nevertheless, considering the crudeness of the approximation, agreement of 

§This is not a test of eq. (80) or eq. (81); this is a test of the quality of the cubic fits. A comparison with the 
direct computation of eq. (81) in Tab. 5 indicates that the coefficients of the linear terms are determined with good 
accuracy, whereas the coefficients of the quadratic terms are determined with about 20% accuracy. 

^Hcre, we mean that the entire Vc+lM is evaluated in the large-/?o approximation, i.e. we set a n = (5/?o/3) n in 
eq. (5) and all f3 n = except (3o in eq. (6). This should not be confused with Vc+l(t) up to NNNLL where (only) 
a,3 is evaluated in the large-/3 approximation. 

l+f. The scale dependence of Vc+l(0 up to NNNLL is by far smaller. 
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the prediction with the lattice data is remarkably good. Vc+l(?") with ats(Q) = 0.2 in Fig. 20 
is much closer to the lattice data than Vc+i,(r) at LL in Fig. 17. On the other hand, a more 
detailed examination reveals limitations of the large-/?o approximation. In Fig. 19, we plot 5£\js( r ) 
computed from Vc+lO*) i n the large-/5 approximation with cts(Q) — 0-2; a cubic fit using the 
data points at rA^° p < 0.5 is shown in the same figure and in Tab. 3. Although the magnitude 
of SE\js(r) is small, 5E\js(r) is oscillatory. The sizes of the coefficients of the polynomial fit are 
unnaturally large, and the fit does not reproduce 5E\js(r) beyond rA^ op = 0.5. Since there 
is no reason to believe that the large-/?o approximation is very close to the infinite-order limit, 
we consider the different behaviors between 5E\js(r) in this case and that of RG analysis (LL 
- NNNLL) to be an indication that the consistency checks performed in this subsection have 
sensitivity to the details. 

5.2 Determinations of roA^g and SE\js(r) in C+L scheme 

From the above analysis, we conclude that all the theoretical expectations derived from OPE are 
either positively confirmed or consistent with lattice data within the present level of uncertainties. 
In this subsection, with the aid of theoretical predictions of OPE, we estimate the infinite-order 
limit of 5E\js(r) in the C+L scheme.* 

Up to this point, we used the central value (0.602) of the relation between r and A^ op in 
eq. (83). Now we examine how our determination of 5Eus(r) will be affected if we vary this relation. 
To simplify the argument, for the moment let us suppose that the lattice data set (when expressed 
in units of ro) has no errors and Vc+l(?") in the infinite-order limit is known. Let V\ at t{f", %) represent 
the lattice data set converted to units of A^i°° p using a given value of x = r A^- op . Then, if x 
equals its precise value Xtrue, according to OPE, 5Eu S (r) = V[ at t(r; x true ) — Vc+lO") goes to zero as 
r —>■ (ignoring an r-independent constant), hence, it would be approximated reasonably well by 

*We note that the analysis in this subsection requires a rather high accuracy in the numerical evaluations of 
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a polynomial of r. When x differs from x trU e, we thus expect 

Vdifr(r; x) = Vi att (r; x) - V c+ h(r) 

= Vi att (r; x tTUC ) - Vc+l(t) + V latt (r; x) - li att (r; x true ) 
« P(r) + A(r i X, Xt rue 

), (84) 

where P{r) is a polynomial of r, and 

A(r;x,x') = V, att (r;x) - V latt {r; x') . (85) 
If we know V[ att (r; x) for some value of x, we can find V\ at t(r; x) for other values of x via 

VW(r;aj')=V htt (r^;x) - r (86) 

On the other hand, we know that VJ att (r;x) tends to const, x (rlog |rAj^g|) _1 at short distances, 
while it tends to a linear potential at large distances. Then, one readily finds that A(r; x, x true ) 
has approximately a "Coulomb" +linear form at rA^ op < 1 if x ^ x truc . 

JV1 o 

In fact, by varying x = r A|^ op within the range given by eq. (83), we find a very good fit of 
V diS (r;x): 

V diS (r;x) ^A(r;x, 0.596) + (0.8 p + 0.7 p 2 ) A^- op . . [NNNLL, a 3 (Pineda)] 

o=rA — ° p 

MS 

(87) 



Inclusion of A(r;x, x') into the fitting function stabilizes the fit considerably when r A^. op is far 
from its optimal value (0.596): The coefficients of the linear and quadratic terms are much less 
affected even if we include cubic term in the fitting function or even if we use data points up to 
larger r for the fit, whereas they are very unstable without A(r;x,x') in the fitting function. A 
further examination shows that it is the inclusion of a Coulombic term (it does not matter very 
much whether log corrections are included or not) that stabilizes the fit when r A^ op is far from 
its optimal value. In principle we should have included A(r;x, x') in the fitting function in the 
previous subsection. A posteriori, it is because the central value of eq. (83) happened to be close 
to the optimal value in eq. (87) that polynomial fits were relatively stable (in particular as we 
increase the order) and all the features seemed quite consistent with the theoretical expectations. 

We note that the error of the input r o^-^° P m eq. (83) is sizable with regard to the accuracy 
of Vc+l(?") in our analysis. Indeed, in [26], the error of the input roA|^ op was the largest source of 
errors in the determination of 5E\js(r). Conversely, this means that our analysis has a sensitivity 
to determine the relation between tq and A^ op by itself (and possibly reduce the error of 5Ens(r) 
at the same time). Hence, in our determination of 5E\js(r), we will determine the value of 
x = r A^ op simultaneously, by performing a fit to the data for ^^(r; x) = Vi att (r; x) — Vc+h{ r )- 

We approximate SE vs (r) by quadratic function (A + A 1 p + A 2 P 2 )Ag op at p = r A|J op < 0.5. 
As for errors, we take into account three types of sources: (i) errors of the lattice data, (ii) error 
of a 3 , and (iii) error due to higher-order corrections. Then we compute the probability density 
distribution for the parameters (Aq,Ai,A 2 ,x) in the following way. Define 



n 

xl = E 

i=i 



VutM; x) - {V c+L (r t ; s) + 1 5V C+L (r t )} - (A + A lPl + A 2 pf )K\ 



op 

MS 



WW 

P v (A ,A l ,A 2 ,x)=Ny 1 J dsdte-^ /2 P s {s)P t {t) ) (89) 
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where x = r A^°° p and = rjA^ op . The normalization constant A/y is chosen such that the 
integral of P V (A , Ai, A 2 , x) over the entire range is unity, s and t parametrize errors of the 
theoretical prediction for Vc+l(V)- Details are as follows. 

(i) We use the first 11 lattice data points given in Tab. 2 of [49]. r\ = rj(x) denotes the distance 
r of the i-th lattice data point (given originally in units of r and r c ) after conversion to 
units of Ag op , using''" x = r Ajg op ; 6Y(x) denotes the error of the i-th lattice data point 

(given originally in units of r ) after conversion to units of A^ op . The first 11 data points 
correspond to pi < 0.5 when x = 0.602. 

(ii) Vc + h(r; s) denotes V c+ h(r) up to NNNLL evaluated with a 3 = s x a 3 (Pineda). We scan s 
between and 2 with equal weight, i.e. P s (s) = l/2if0<s<2 and P s (s) = otherwise. 
The interval < s < 2 covers within its range the estimates of 03 by Pineda [26], by 
Chishtie-Elias [13] and by large-/5o approximation. 

(iii) tSVc+hir) represents an estimate of the difference between Vc+l(?") up to NNNLL and 
Vc+h(f) in the infinte-order limit. SVc+l^) is estimated by the difference between Vc+l(?") 
up to NNLL and that up to NNNLL. We scan t between —1 and 1 with equal weight, i.e. 
P t (t) = 1/2 if |*| < 1 and P t (t) = otherwise. 

Alternatively we may use the QCD force -Fqcd(^) = dVQcr>/dr in the determination of 5E\js(r) 
and x. Since we are not interested in the r-independent part of the potential, we can extract 
information on the relevant parameters using the force as well. Similarly to before, we define 
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i=i 



'F laU ( ri ;x) - {V£ +L (r i ;s)+t5V£ +L (r i )} - (A 1 + 2A 2Pl )A^ 



(90) 



P F (A U A 2 , x) = Afp 1 J ds dt e~ x ? /2 P s (s) P t (t). (91) 

We use the first 10 points of the lattice data for {Fi att (ri), Sf (n)}, corresponding to pi(0.602) < 0.5, 
given in the same table (Tab. 2) in [49]. Other details are same as in the case using the potential. 

There is a considerable difference between the use of the potential and the force in the deter- 
mination of 5Eus(r) and x. The difference stems from different correlations of the errors ({<^}, 
of the respective lattice data sets and from our treatment of these errors. It is known that 
there exists a high correlation among the errors of the lattice data at different r^, for the QCD 
potential or for the force. It is also known that the error correlation of the force is smaller than 
that of the potential [51]. On the other hand, up to now, the covariance matrix of the errors for 
neither of these quantities is available. We therefore decided to use the lattice data with Gaussian 
errors, neglecting the correlations; see eqs. (88) and (90). This treatment should, in general, result 
in overestimates of errors in the determination of (Ai, A 2 ,x). 

Bounds on 5Eus(r) can be obtained from the probability density distributions for (Ai,A 2 ), 
defined by 

Pl lM {A u A 2 ) = J dxdA P v (A ,A l7 A 2 ,x), (92) 
Pl M {A u A 2 ) = [dx P F (A u A 2 ,x). (93) 



tWe fix r c /r = 0.5133 [49]; its error is small, which we neglect in our analysis. 
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Using the QCD Potential 
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Figure 21: (i) Contour plot of the probability density distribution » (A\, A2), corresponding to 
68% and 95% CL regions. Cross represents (Ai,A2) with the highest probability density, (A, A) = 
(0.40,0.76). (ii) Bounds on 5E\js(r) in the C+L scheme corresponding to the regions of (i). Quadratic 
fit with the highest probability density is also plotted (blue line). 



Using the QCD Force 
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Figure 22: (i) Contour plot of the probability density distribution Pj[ A (A±, A2), corresponding to 
68% and 95% CL regions. Cross represents (A, A) with the highest probability density, (A, A) = 
(1.03,0.60). (ii) Bounds on 5E\js(r) in the C+L scheme corresponding to the regions of (i). Quadratic 
fit with the highest probability density is also plotted (blue line). 



(ii) 
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Figs. 21,22(i) show contour plots of these probability density distributions corresponding to the 
68% and 95% confidence level (CL) regions. The corresponding bounds on 5i?us( r ) i n the C+L 
scheme are given by [see Figs. 21, 22(h)] 



Using the potential: 

-0.7p+ l.Op 2 < 5E vs (r)/A?± op < 2.1p + 0.2p 2 (68% CL) 



MS 

-1.3 p + 1.1 p 2 < 5£ us (r)/A|± oop < 3.7p-0.8p 2 (95% CL) 



(94) 



Using the force: 

0.2p + 0.7p 2 < 5£ us (r)/A|J op < 2.0p + 0.4p 2 (68% CL) 

(95) 

-0.2p + 0.8p 2 < 5£ us (r)/A|J op < 2.6p + 0.2p 2 (95% CL) 



These bounds are mutually consistent, and the latter bounds are tighter. Since the origin 
(Ai,A 2 ) = (0,0) lies outside the 95% CL regions in Figs. 21,22(i), we conclude that 5Eus(r) 
being 0(AQ CD r 3 ) or 5Eu S (r) = is disfavored. We see that positive A 2 is favored, in agreement 
with the fits in Tab. 3. 

The probability density distributions for x = r A|^ op are defined by 

Px(x) = J dA dA 1 dA 2 P V (A , A u A 2 ,x), (96) 
Pf (x) = J dA x dA 2 P F (A l ,A 2 ,x). (97) 

They are shown in Fig. 23. Each distribution is close to Gaussian, so we simply quote the mean 
and the standard deviation for x: 

Using the potential: x = 0.592 ± 0.062 , (98) 

Using the force: x = 0.574 ± 0.042 . (99) 
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Note that we did not use the relation (83) at all to obtain these results. Our results are mutually 
consistent, as well as in excellent agreement with eq. (83). * The error in eq. (99) is of similar size 
to (slightly smaller than) that in eq. (83). 
Some comments are in order. 

As explained above, the sensitivity to x = A^ op originates from the mixing of a Coulombic 
term in V<nff(r;x) when x is different from its true value. The only assumption we made in 
the determination of x is that 5Eys( r ) can be approximated by a quadratic polynomial of r at 
r A|^ op < 0.5. Since 5E\j S (r) goes to zero at sufficiently small r according to OPE, a polynomial 
fit should be reasonable. Effects of higher powers of r is expected to be suppressed at small 
r. In fact, we made a consistency check by including A 3 p 3 term into the fits of SEysi^) and/or 
by varying the number of data points used for the fits. As we include more data points, errors 
determined from P^{x) and -Pj(x) decrease, respectively, while the quality of the quadratic fits 
tends to get worse (quality is better with the cubic fits). In any case, the obtained bounds on x 
are consistent with eqs. (98) and (99). We found that our present choice, quadratic fit with the 
first 11 (10) data points, is close to optimal in performing the fits. 

When we scan s and t between the intervals (0, 2) and (—1, 1), respectively, the minimum values 
of Xv an d Xf var Y between (0.5, 0.6) and (0.1, 0.2), respectively. Since the number of degrees of 
freedom is N doi = 11 - 4 = 10 - 3 = 7, both (xv)mm/N do{ and (x|0min/N d of are below 10%. This 
is consistent with existence of high correlations among the lattice errors at different rj, which we 
mentioned already. We note that our errors in eqs. (94), (95), (98), (99) may be overestimated for 
this reason. 

The difference between the bounds obtained by using the potential and force can be attributed 
to the difference of the correlations of the lattice errors. Since the correlation is larger for the 
potential, the errors of the lattice data are effectively more enhanced (overestimated) in our 
treatment, hence the bounds are wider when we use the potential. The lattice errors are the 
dominant source of errors in the determination of (A\, A2) and x, both when we use the potential 
and force. In this sense, the covariance matrices of the lattice data are highly demanded. 

Since the errors for Vc+h( r ) (parametrized by s and t) are much larger than the errors of the 
lattice data, one may wonder why the latter can be dominant source of errors with regard to the 
former. [Note that in Fig. 19 the errors of the lattice data are smaller or comparable to the size of 
the symbols used for the plot; the variation of Vc+l(^; s) with s is comparable in size to SVc+h(r).] 
This can be understood as follows: 

(1) Practically, the measurement of x is sensitive only to the "Coulomb" part of V d is(r;x), hence 
only the "Coulomb" part of the errors matters. Let us denote by V^ L (r) the difference between 
Vc+lM up to N n LL and V c+h (r) up to N n_1 LL. Then we perform fits of V^ L (r)/A^ op in the 
form c-i p^ 1 + C1P + C2 p 2 , using the 11 data points evaluated at r = r-j. The results for n — 1, 2, 3 
are shown in Tab. 4. Magnitudes of all the coefficients decrease as the order increases, if we 
take s = 1 as a reference for V^) h (r). This is natural, since Vc+ L (r) — > as n — > 00, therefore all 
— * 0. Thus, it is quite reasonable to estimate the "Coulomb" part of the errors of Vc+l(?") by 
the "Coulomb" part included in t V^ L (r; s = 1) or by that included in V^ L (r; s) — V^ L (r; s = 1). 
In fact, these error estimates are encoded in our analysis. Noting that we neglect the correlation 
of the errors of the lattice data, the errors of the lattice data are indeed larger than the "Coulomb" 
part of the errors of Vc+l(?")- [Of course, the argument given here is only for demonstration to 
understand the errors better. When we derived our results eqs. (94), (95), (98), (99), neither did we 
do a fit of the form c_i + c\ p + C2 p 2 nor extract a "Coulomb" part.] 

*We also note the value obtained by [49], 0.586 ± 0.048. This value is closer to our values. 
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KSlM/A^ -O.O^Sp^ + l.lp-O.Sp 2 

^ 2 + } L (r)/A|± oop -0.0087 p- 1 + 0.8p - 0.2 p 2 

l/g L ( r; s = 1)/A^° P -0.0028 p' 1 + 0.5 p- 0.05 p 2 

[V^ L (r; 5 = 2) - V^Jr; s = l)]/A|± oop -0.0038 + 0.6 p - 0.08 p 2 
[v£ } L (r; 8 = 0) - g = 1)]/AS° P +0.0038 p' 1 - 0.6 p + 0.08 p 2 

Table 4: Fits of Vc™L(rj)/A^ op m the form c_i p _1 +Ci p+c 2 p 2 , using the 11 data points evaluated 
at r = 7^(0.602). 



(2) There is a similar mechanism in the determination of (Ai,A 2 ). Since A(r;x,x') has a 
"Coulomb" +linear form, if a small admixture of "Coulomb" part is allowed in Vdiff(r) due to 
the errors of the lattice data, the linear term attached to the "Coulomb" part mixes in as well. 
Hence, if the lattice errors are larger, allowing a "Coulomb" part to mix in, the bounds on (Ai, A 2 ) 
spread mainly in the A\ direction. This explains the difference of Figs. 21 (i) and 22(i). Again, it 
is the size of the "Coulomb" part of the errors that matters. 

Our results can be compared with the determination of 5Eu S (r) by Pineda [26]. It is the only 
study that determined the non-perturbative contribution using OPE, preceding our current work. 
There are some important differences between Pineda's analysis and ours. 

• Pineda used x = r A^ op as an input parameter, given by eq. (83). Its error turns out to 
be the dominant source of errors in the determination of 5E\js(r) (defined in RS scheme 
[14, 26]). On the other hand, in our analysis, we determine x from a fit to the data. 

• We estimate the error of the singlet potential (in the C+L scheme) by varying s and t 
between < s < 2 and \t\ < 1. On the other hand, Pineda estimates the error of the singlet 
potential (in RS scheme) by varying 03 between the range corresponding to 1/2 < s < 3/2, 
while there is no estimate corresponding to variation of t, i.e. t is fixed to zero. Thus, our 
error estimate of the singlet potential is more conservative. 

• Pineda does not incorporate errors of the lattice data at all. On the other hand, in our 
analysis, they are included neglecting the correlation. Since our analysis is sensitive to a 
"Coulomb" part, the lattice errors are the major source of errors. 

• The singlet potential in RS scheme contains (9(AQ CD r 2 ) renormalon. It means that the sin- 
glet potential has an intrinsic uncertainty of this order. This essentially prevents a determi- 
nation of 5Eu S (r) with better than 0(AQ CD r 2 ) accuracy, because 5Eu S (r) = Vqcd(^) — Vs(r) 
cannot be defined with better accuracy. On the other hand, our potential is free from the 
(9(AQ CD r 2 ) renormalon (and also from the rest of IR renormalons) . So, at least conceptually, 
there is a difference in the achievable accuracies between the two analyses. 

Due to these differences, comparisons between our bounds on 5E\js(r) and those of Pineda are not 
straightforward. An explicit bound obtained by Pineda, assuming a form 5Eus(r) = const, x r 2 , 
reads 

\SE vs {r)\ < 2.8 (A^ op ) 3 r 2 (RS scheme, Pineda [26]). (100) 
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[Vc +L (r) - 




LL 


0.53 p 2 


+ 0.00 p 3 + • • • 


NLL 


0.85 p 2 


+ 0.07p 3 + ••• 


NNLL 


1.04 p 2 


+ 0.27p 3 + ••• 


NNNLL 


1.02 p 2 


+ 0.58p 3 + ••• 



Table 5: Expansion of [Vq+lM - V^ R) (r; M/ )]/Ag op in p = rAg op , computed using eq. (81). We 
neglect p-independent constants. V^ R ^(r;/ij) is computed in the second scheme and with p.f = 3Aj^° p . 
Other parameters for Vc+l(0 and Vg R \r;pf) are same as in Fig. 17. 

One way of comparison may be to perform a fit by setting A\ = in our analysis. From the 
probability density distribution -Pj^/^O, A 2 ), we obtain 

n n / a 3-1oo P n3 2^x7-1 / \ ^ n n /a3-1oo P \3 2 ( C+L Scheme, Ai = \ /im\ 

0.7 (A™ p ) r < d_E us (r) < 0.9 (A™. p ) r _ o0y „ T . . . (101) 

v ms ^ ^ ms ^ ^ 68% CL using force / v ' 

5.3 Determination of ^E^usC^) in factorization scheme 

We can also determine the size of SEu S (r; pf) in the factorization scheme. In fact, in the second 
scheme (within the factorization scheme), the purely non-perturbative contribution is common to 
that in the C+L scheme. This is because the difference of 5E\js(r) is merely the difference of the 
Wilson coefficients Vc+h(f) and V^^r; pf); it is given by (minus) the right-hand-side of eq. (81), 
which is systematically computable by means of perturbative expansion and log resummation via 
RG. So, our task reduces to estimating the infinit-order limit of eq. (81). Using eq. (81), we 
calculate Taylor expansions of Vc + L(r) — Vg K \r;pf) in r, which are listed in Tab. 5 for pf = 
3 A^ op . We estimate errors by the size of the NNNLL corrections and obtain 

MS J 

[V c+h (r) - Vf } (r; 3 A|L oop )]/Ag op » (1.02 ± 0.02) p 2 + (0.58 ± 0.31) p\ (p < 0.5) (102) 

Thus, 5Eu S (r; 3 A|^ op ) in the factorization scheme (second scheme) is given by the sum of eqs. (95) 

and (102). Given the above estimate, estimating Vc+l(?") — Vg (r;pf) for other value of p/ is 
straightforward, since the difference of the integrals [eq. (81)] can be evaluated by integrating over 
q along the real axis and the integrand is free from singularities; convergence is fairly good in this 
region of q, so one may simply use the prediction up to NNNLL. 

By the same token, we can estimate 5E\js(r) in the first scheme within the factorization 
scheme. The difference between the first scheme and second scheme is perturbatively computable. 
In practice, up to NNNLL, we do not find a significant difference between the first scheme and the 
second scheme. We consider that we may apply the above estimate (102) also to the first scheme 
within the factorization scheme. 

6 Summary and Conclusions 

In this paper, we analyzed the static QCD potential in the distance region relevant to heavy 
quarkonium spectroscopy, 0.5 GeV _1 (0.1 fm) < r < 5 GeV _1 (l fm), using perturbative expansion 
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and OPE as basic theoretical tools. The analysis consists of three major steps: 

(I) Behavior of the QCD potential at large orders of perturbative expansion was analyzed. As 
for the higher-order terms, we used the estimates by large-/3o approximation or by RG equation; 
as for the renormalization scale /i, we varied it around the minimal-sensitivity scale [or, more 
precisely the scale defined by eq. (21)]. Then the perturbative expansion of the QCD potential, 
truncated at 0(a$ ), was separated into a scale-independent (prescription-independent) part and 
scale-dependent (prescription-dependent) part when N >> 1: 

V N (r) = V c {r) + B(N, +Cr + V(r, N, £) + (terms that vanish as N -> oo). (103) 

Here, £ is a parameter for changing the scale (£ = 1 corresponds to an optimal choice of scale). 
Vc{r) is a "Coulomb" potential, which includes logarithmic corrections at short-distances; B is 
an r-independent constant; Cr is a linear potential; T>{r) behaves as r 3 ^ -1 x (log corr.) if £ < 1, 
whereas it is 0(r 2 ) and divergent as N — > oo if £ > 1. Vc(r) and Cr correspond to a renormalon- 
free part of V/v(r) and are finite and independent of £; thus the scale-independent part has a 
"Coulomb" +linear form. On the other hand, B and V(r) correspond, respectively, to the (9(Aq C d) 
IR renormalon and beyond (9(Aqcd) IR renormalons (starting from the (9(AQ CD r 2 ) renormalon) 
contained in V/v(r); they are dependent on £ and divergent as N — > oo if £ is sufficiently large. 
Detailed analytic behaviors of each component have been studied. 

(II) In the framework of OPE of the QCD potential, (a) we gave explicit renormalization pre- 
scriptions for the Wilson coefficient (singlet potential) Vg {r;/if), which belong to the class of 
conventional factorization schemes with a hard cutoff; (b) the scale-independent part of (I) was 
generalized and promoted to a Wilson coefficient Vc+l(?"), which is independent of 'the factorization 
scale Hf. Both Vg R \r; and V c +l(^) are free from IR renormalons and IR divergences. Several 
properties of these Wilson coefficients and of the corresponding non-perturbative contributions 
have been derived (partly already in [24]): 

• ^ / 5 R ' ) ( r !/ i /) an d Vc+hir) can be computed systematically using perturbative expansion and 
log resummation via RG, and (in principle) the predictions can be improved to arbitrary pre- 
cision. Hence, the corresponding non-perturbative contributions 5Eu S (r) are unambiguously 
defined. 

• With the usual hierarchy condition Aqcd <C A*/ ^ l/ r ) the difference between Vg R \r;/j,f) 
and Vc+h( r ) is 0(fi 3 jr 2 ) and perturbatively computable. Vc+l(?") is closer to VqcdM than 

• 5E\js(r) in the factorization scheme is 0(fi^r 3 ) at very short-distances (AV(r) ^> ///), 
whereas it is 0{ix^r 2 ) in a semi-short-distance region (AV(r) <C /// 1/r). 

• 5i?us( r ) i n the C+L scheme (/^/-independent scheme) is (9(AQ CD r 3 ) at very short-distances 
(AV(r) ^> Aqcd), whereas its behavior cannot be predicted model-independently in a semi- 
short-distance region (AV(r) ~ Aqcd)- 

We conjectured that the region of our interest corresponds to a semi-short-distance region where 
AV(r) ~ A QCD (<A*/)- 

(III) We computed V^ R ^(r; /x/) and Vc+l(^) numerically for n; = at our current best knowledge 
(NNNLL with certain estimates of a 3 )* and compared them with the lattice computations of 

*As far as US logs are concerned, we observe that their effects are small, hence, we did not resum the US logs 
but included only up to NNNLO in the analysis. 
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Vqcd(^) m the quenched approximation. We confirmed that the theoretical predictions of (II) are 
either correct or consistent within the present level of uncertainties. * We find that a linear potential 
in 5Eus(r) reduces with increasing order, consistently with vanishing in the infinite-order limit; 
at NNNLL, it is much smaller than the string tension as determined by lattice simulations. Then, 
we performed fits of Vdm{r; x) = VJ a tt(r; x) — Vc+l(?") and determined simultaneously 5E\js(r) and 
x = r A^ op (relation between Sommer scale and A^g). A sensitivity to x originates from mixing 
of a Coulombic term into Vdis(r;x) when x differs from its true value. Both the QCD potential 
and QCD force were used for the fits. The latter resulted in tighter bounds due to a smaller 
correlation of lattice errors. We obtained 

r Ag op = 0.574 ± 0.042 , (104) 

in excellent agreement with the determination via Schrodinger functional method, eq. (83) [47]. 
We also obtained 

0.2p + 0.7p 2 < 5E vs (r)/A^ op < 2.0p + 0.4p 2 , (C+L scheme) (105) 

where p = rA^ op . [See also bounds on the coefficients of quadratic polynomial in Figs. 21,22(i).] 
In the factorization scheme, we obtain, for instance, 

no li 7 2/rp / WA 3-Ioop/ On , , A 2 ( factorization scheme \ 
0.2p+1.7p < 5E V s{r;fi f )/A— < 2.0p + 1.4 p . ^.s-ioop (106) 

Estimating SE\js(r; pf) for other pf is easy. In the factorization scheme, the obtained bounds 
are consistent with 0(p 3 jr 2 ), rather than 0(pjr 3 ). In the C+L scheme, the obtained bound 
is consistent with (9(AQ CD r 2 ) (vanishing linear potential) at 95% CL, but existence of a small 
linear term is more favored; furthermore, 5E\js(r) ~ (9(AQ CD r 3 ) or 5E\js(r) = is disfavored. 
Consequently, we find that pj ^> AV(r), in accord with our conjecture. Also AV(r) ~ A QCD is 
more likely than AV(r) 3> Aqcd.* 

The analysis (I) provides a reasoning within perturbative QCD, why we observed agreement 
between the recent perturbative computations of the QCD potential with phenomenological po- 
tentials or lattice results: in the large-order limit, Vjy(r) does approach a "Coulomb" +linear form. 
It is quite intriguing that we can separate the renormalon-free part and renormalon-dominant 
part in a natural way. Furthermore, in this analysis, the coefficient of the linear potential can be 
computed analytically up to NNLL. 

In (II) we defined renormalization prescriptions for Vg (r;pf), in which all the known IR 
renormalons are subtracted. § Moreover, introduction of a factorization-scale independent scheme 
for the Wilson coefficient, Vc+h(r), is new. Vc+l(?") has some appealing theoretical features: it has 
no intrinsic uncertainties, is systematically computable, and is closer to Vqcd(^) than Vg R \r; pf). 
In any case, since the differences between the different schemes are perturbatively computable 
with good accuracy, the C+L scheme would serve as a useful reference. 

^We also observed some limitations of the large-/?o approximation. 

*Although we have not considered the possibility AV(r) <C Aqcd in this paper (since it seems to lie outside 
the applicable range of our analysis), it may be worth examining this possibility in detail in view of our present 
results. 

§We ignored the instanton-induced renormalon singularities, which are known to give very small contributions. 
From a general argument, there should also exist contributions to IR renormalons, which cannot be written in the 
form of eq. (18); we neglected them too. 
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In the numerical analysis (III), the only assumption we made in the fits is that we can ap- 
proximate 5E us (r) by a quadratic polynomial of r at rA|^ op < 0.5. This should be reasonable 
since 8Ejjs(r) — > as r — > according to OPE. (An r-independent constant is irrelevant in our 
analysis.) We checked validity of the assumption by including cubic term and by varying the 
number of data points used for the fits. 

In Fig. 17 we have seen apparent convergence of Vc+i,(r) towards the lattice data, up to our 
current best knowledge. However, a closer examination of 8E\j S (r) in the C+L scheme revealed 
that SE\js(r) = is disfavored. In fact, we have improved both quantitatively and conceptually (in 
the sense that we removed all the renormalons from the singlet potential) the bounds on 8Ejjs(r), 
as compared to the pioneering study by Pineda [26]. 

The OPE analysis of the QCD potential provided, as a byproduct, a new method for deter- 
mining x = r A^ op . Our current result gives an error comparable in size to the error of the 
conventional result using the Schrodinger functional method [47]. The mechanism for the sensi- 
tivity is fairly clear, as well as sources of errors are understood well. The present status is that 
the errors of the lattice data contribute more significantly than the errors of Vc+l(?")- Hence, 
information on the correlation of the lattice errors (in particular, the covariant matrix) is highly 
demanded in order to reduce the error. 

Finally let us comment on the applicable range of perturbative expansion and OPE of VqcdO")- 
We saw in Fig. 17 that the current best perturbative prediction of the Wilson coefficient Vc+l(?") 
follows the lattice data up to r < r ~ 0.5 fm. We consider that the distance, at which string 
breaking occurs, serves as a measure of the distance where the perturbative expansion breaks 
down (in the theory with rii > 0). It is around 1 fm according to the recent lattice simulation 
[52]. It is clear that the string breaking phenomenon is non-perturbative and that the present 
perturbative computation of the QCD potential lacks ingredients necessary for description of this 
phenomenon. In the context of heavy quarkonium phenomenology, string breaking corresponds to 
the decay T(45 l ) — > BB. Since empirically the root-mean-square radius of T(4S) is around 1 fm, it 
is consistent with the lattice results. We also know empirically that phenomenological potentials 
are approximated well by a Coulomb+linear form at r < 1 fm. It means that, if we separate 
heavy quark and antiquark, we have a sensitivity to the linear potential at distances before string 
breaking takes place. Thus, we consider r < 1 fm [corresponding to heavy quarkonium states 
below T(4S*)] to be the range in which perturbative expansion may make sense. (Certainly more 
terms of the perturbative expansion need to be included in order to have an accurate prediction 
of Vc+l(?") as r approaches 1 fm.) Already in this range the QCD potential exhibits a linear 
behavior in addition to the "Coulomb" part. Our analysis indicates that the qualitative argument 
presented at the end of Sec. 2.4 may be valid in this very range. Ultimately, it depends on whether 
the linear term in 5E\js(r) is truly vanishing or not." 
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^Note that the linear potential we are concerned here has, a priori, nothing to do with the linear potential at 
r >• AqJ, d usually associated with confinement (in the theory with ni = 0). A priori, we see no reason that both 
linear potentials should have a common slope. Nevertheless, empirically these two linear potentials seem to have 
a common slope. 
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Appendices 



A Basic Formulas 

In this Appendix we present detailed formulas useful for computing a v T (q) up to NNNLO as well 
as V)v(r) for a finite but large N. 

The perturbative expansion of the l^-scheme coupling in momentum space a v T (q) is defined 
by eqs. (2)-(5). The polynomials in eq. (4) up to NNNLO are given by 

Po(£) = ao, (107) 

P 1 (£)=a 1 + 2a o o £, (108) 

P2(£) =a 2 + (Aa 1 o + 2a o 1 )£ + 4a o o 2 £ 2 , (109) 
Ps(t) = % + (6 a 2 p + 4 ai A + 2 a &) ^ + (12 a x A, 2 + 10 a ft A) £ 2 + 8 a ft 3 f 5 , (110) 

where 

t = \og{n/q). (Ill) 
The coefficients of the beta function (3 n , defined by eq. (6), are given explicitly by 

ft = ll-^n,, (112) 
/3 1 = 102-yn,, (113) 
A = !fl - 5™, + f„ ? [53,54,, (114) 

A =2^5? + 3564 6 
6 

/ 1078361 6508 ( 3 \ / 50065 6472&\ 2 109JW 

+ 162 27- J n ' + (, "W + "81- J n ' + -729- [55] ' (U5) 

where (3 = C(3) = 1.2020... denotes the Riemann zeta function ((z) = J2^=i V 71 * evaluated at 
z = 3. 

Presently, known up to n = 2. 

o = l, (116) 
31 10 , , 

°i = -3 " Y Hl [1 ' 2] ' (U7) 

4343 ^ 2^ RR , 9vr 4 /1229 52C 3 \ 100 2 
a 2 = ^- + 36tt +66(3-^- ( ^ + ^- I + — ^ [4,5]. (118) 

It is known that 03 is IR divergent; the coefficient of the divergence and associated logarithm have 
been computed [10, 11]: 



a 3 = 72tt 2 



^ + 4{2£ + log(47r)- 7s } 



as, (H9) 



where the IR divergence is regularized by dimensional regularization (D = 4 — 2e). a 3 is just 
a constant independent of e, /i, r. So far, only some estimates of its size are known: e.g. 
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a 3 (large-A)) = (|/3 ) 3 ~ 6162, a 3 (Pineda) w 18688 [14, 26], a 3 (Chishtie-Elias) « 20032 [13] 
for n/ = 0. 

Physical logarithm associated with the IR divergence can be extracted as follows [10, 24]. 
Instead of defining the ultrasoft contribution 5Eus(r) m ec i- (17) as a non-perturbative quantity, 
it can be computed in a double expansion in as and logas within pNRQCD:* 

CVC> 5 (/i) 4 



6 Eh 



us 



(r)} 



double 
cxp. 



247rr 



X 



X - + 81og(/ir) - 21og(o s (//)) + | + 2 lE + 41og(4vr) 



+ o(4), 

(120) 

where 7# = 0.5772... denotes the Euler constant. Upon Fourier transform, 1/e and log// terms 
of eqs. (119) and (120) cancel each other. 1 ' The remaining log(CAO; s ) is the physical logarithm. 
The argument of the logarithm in eq. (120), Cac^s — 2rAV(r), represents the ratio of the IR 
regulator 2AV(r) and 1/r; cf. Sec. 2.4. If we perform OPE in conventional factorization schemes, 
we introduce the factorization scale /if, and the IR regulator 2AU(r) will be replaced by [if. Thus, 
one finds the ultrasoft logarithm 

= _ c F e> 5 (^ 

us-iog 24 7rr 6VP/ ; v ; 

It is easy to verify that the as(fi) 4 log(/i//g) term of eq. (7) generates eq. (121) after Fourier 
transform. (Oppositely, using RG equation with respect to evolution of jif within pNRQCD, one 
can resum US logs as given in eq. (7) [12].) 

One may verify explicitly that a v T (q) up to NNNLO, defined via eqs. (107)-(119), and 
[SEus(r)] de in eq. (120) are separately consistent with RG equations with respect to evolution of 



/^ + /^))^y 



X = 



X = a v '(q) or 5E m (r) 



double 
cxp. 



(122) 



This should be so, as long as Vqcd(^) and cty T (q) are defined from the Wilson loop via eqs. (2)-(4), 
since the Wilson loop is independent of fi} To verify eq. (122), one should note that the beta 
function in general dimension (in MS scheme) has a form 

= -tots + &?(as)Lo • (123) 

In the rest of this Appendix, we present formulas useful for computing Vjv(r) for a large (but 
finite) N. The expansion of as(q) n in terms of cxs(fj) can be obtained by iterative operation of a 
derivative operator as 



®s(q) r 



exp 



E 

fc=0 



X 



-2£) 



k r 



k\ 



(3{x) 



d_ 

dx 



i k 



X n 



(124) 



*It can also be obtained from the difference between the resummation of diagrams in Fig. 1 and its expansion 
in as before loop integration. 

'''Note that the expansion of VqcdM m as, obtained from a v T (q), coincides with the expansion of the bare 
singlet potential Vs(r) in as; cf. Sec. 2.4. Hence, the sum of [<5.Eus (?*)] d c , as given by eq. (120), and Vs(r) 
represents the expansion of Vqct>(t) in as and log as- 

tThere exists a definition of the singlet potential through threshold expansion of diagrams contributing to a 
quark-antiquark Green function [56]; in this definition, ^-independence of the potential is not preserved. 
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C s 



L 8 



Figure 24: Integral contour C s . 

where f3(cxs) is defined in eq. (6) [and eq. (123) if necessary]. 

The V-scheme coupling in position space, ay(l/r), is defined by Vqcd(?") = _ Cf ay(l/V)/r. 
The series expansion of «y(l/r) in terms of the MS coupling renormalized at \x = exp(— 7_e)/V is 
obtained as follows [57]. Using the coefficients g m defined by 



m=0 



m u m = exp 



CO u / , s J. 

£«fil {2 ..l.(. 1) . } 



fc=2 



(125) 



we may write 



«v(l/0 = J^#r, 



m=0 



-(3{x) 



d_ 

dx 



n=0 



a, 



(4tt)' 



(126) 



x^a s (exp(-~f E )/r) 



To obtain the expansion of dy(l/r) in terms of as(/x), we first compute eq. (126), then substitute 
the expansions of as(exp(— 7£)/r) n in terms of cxs(fj) computed using eq. (124). By truncating 
the series at an appropriate order in as(/x), the result reduces to a polynomial of log(/ir). 



B Integral Representation of [cts(q)]oo at NLL 

[as (q) ]oo for the 2- loop running coupling constant can be expressed in a one-parameter integral 
form. After integrating the RG equation, as(q) is given implicitly by the relation 



log (o/A— op , 



2tt 5 / 4tt . . 
f - log ( - r <) ) 



Joas 



(127) 



Hence, using Cauchy's theorem, one may write 



MQ) = jj c ds (s)' 1 [log (?/Aj£*) + S - log (-2es) + s 

- £ ds (sy 1 J™ dx exp -x | log (g/Ag° P ) + | log (-2es) + s 



(128) 



A 



o Jo 



* (<f/ A S° P ) 



I\^/2 1 

2iJ r(l + x<J/2) 



(129) 



The integral contour C s is shown in Fig. 24. In the last equality, we rescaled s — > s/x and used 
the formula l/r(z) = i(2vr)~ 1 / c ds {-s)~ z e~ s . 
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The truncated series expansion of as(q) in a M = a,s{fi) can be obtained as follows. One rewrites 



A^ op in terms of in eq. (128). After changing variables as s 



' 47T 



+ 5)u, we have 



du (—u) 1 



1 + 



u 



a^S + 47r//3 ' 27r 



+ 



log(^) + I log(-eu) 



(130) 



Expanding the integrand in and integrating at each order of the expansion, one obtains the 
series expansion of as(q) in a^. It is then straightforward to truncate at order a^ N . Sending 
N — > oo, we obtain 



lim 



a. 



1 + 



u 



+ 



a^ + An/po ' 2tt 
5 



l0g (f ) + 2 log(_eM ' 



A? 



l0 S (^/ A S° P ) + 2 log ^ 2eS ) + S 



x 1 - exp 



-3e{log(g/A£?° P )+^log(-2e S ) + S 



(131) 



where we re-expressed the truncated series in terms of s and A^ op . Similarly to eq. (129), we 
find 



2vr 



3* 



[asiq)U = fol dX ( g/A S° P ) X ^) XS/2 T(l + x5/2y 



(132) 



Thus, the one-parameter integral forms given in eq. (52) are obtained. 



C Analytic Formula of the Linear Potential in Case (c) 

We present the analytic formula for the coefficient of the linear potential, defined by eq. (39), 
in case (c). The integral eq. (39) can be reduced to a one-parameter integral form by change of 
variables, e.g. from q to z = 1/cxs- Then one readily sees that the integral can be expressed in 
terms of the confluent hypergeometric function except for the coefficient of ct2, while the coefficient 
of a 2 can be expressed in terms of generalized confluent hypergeometic functions. 

For convenience, we first define some auxiliary parameters. The two solutions to the quadratic 
equation 



a s 



5>(2) n=0 (133) 



\ 47T 

case(c) n= Q 



[cf. eq. (6)] are denoted as§ as = and 1/uj*. Then we define 



(134) 

p OJ — UJ* 



We also write bo = /3q/(Att). 



§ We assume that the two solutions are complex conjugate of each other. This is the case when the number of 
active quark flavors is less than 6. 
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The coefficient of the linear potential in case (c) is given by 



where 



✓7(c) _ 2ttC F f A 3-loop\ 

3 n v a hs ; 



2 



a n + ^-n 1 + -^-r n 2 

4:71 (47Tj 2 



(135) 



n = 2Re[F(p+ 1 1 ,p*)+LuF(p,p*)] , (136) 

ll 1 = 2Re[F(p,p*)}, (137) 

TZ 2 = 2Re[G] -b s {-u) 2p - 1 {-u*) 2p '- 1 . (138) 

The function F is defined by 

F(x,y) = br y+&K - J r(1 _ 2x) W y _ x , x+y _ 1/2 (^), (139) 

in terms of the Whittaker function, which is related to the confluent hypergeo metric function ±F\ 

as 

W^{z) = r (~ 2 ^ ^+V2 e -,/2 iFi{pL _ « + l, 2„ + 1; s) 

1 ll ~~ z 1 ~ K ) 

I ( 2/- ) _ M+1/2 e _, /2 |7 .. ( ;/ _ K+ i_ 2/Jl+ 1; ^ (IK,) 



r(| + / .-«) 

On the other hand, G is defined by 

G = b 5 e («V&o)+**(2p'-l) 



x I 7 t~ r r\ f 1, 2p - 1, 2 + <S; ^) 

1(1- 2p) 5(1 - 2p, 1 - 2p*) V ' ^ ' ' 6 « / 

l -2-5 
00 



— sin(27rp) r(-2 - 5) ~ 2 (l, 1 - 2p*, 3 + 5; ^) }. (141) 
Ti and E 2 represent Appell confluent hypergeometric functions [58] defined by the double series 

< J ml n 1 



ml nl 

m,n=0 

00 



^•™»>=E<#7&*v- < 143 > 

where (a) n = T(a + n)/r(a) is the Pochhammer symbol. 

D Numerical Evaluation of V^ R ^(r;/xj) and Vc+lO*) 

In this Appendix, we present a method for accurate numerical evaluations of Vg R \r;/j,f) and 
Vc+l(?")- The former is defined in Sec. 4.1 to be 

t/( r V ^ 2C F [°° sin(gr) (R) 

( r ;/ i /) = / d( l Uvs WiHf) (144) 

J (if Q r 
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with 



avJiWf) = a v T ((l) + Sa^q; li j) (145) 

N , s 



n=0 



N = 0,1, 2, and 3 correspond to V^ R) (r; /i/) up to LL, NLL, NNLL and NNNLL, respectively. (We 
do not resum US logs but include only up to NNNLO.) = P n (0) = a n for n < 2, whose explicit 
forms are given in eqs. (116)— (118); a^ s = a 3 + 144-7T 2 log(/i//g) in the first scheme [corresponding 
to eq. (74)], while a^ s = a 3 — 1207r 2 + 1447r 2 {7^ + log[3as(g)]} in the second scheme [corresponding 
to eq. (75)]. 

We deform the integral path of eq. (144) into upper half plane: 



r\/~i POD 

V^ R \r; p f ) = ^ Im / dk 

n Jo 



exp(iqr) (R) ' 



qr 



(147) 



In order to evaluate the integral numerically, we first solve the RG equation (6) with a given input 
value (e.g. ets(Q) — 0.2) and find the value of Q/A^ and the value of as(fJ>f) for a given ////Ajjjg. 
Then we solve the RG equation (6) along the integral path q — /if + ik (0 < k < 00) in the 
complex plane, with ag(fif) as the initial value. In solving the RG equation, we take the sum for 

n < 0, 1, 2 and 3 on the right- hand- side of eq. (6), respectively, corresponding to Vg K \r; ///) up to 
LL, NLL, NNLL and NNNLL. 

The singlet potential in the /^-independent scheme Vc+i,(r) is defined in Sees. 4.1 and 4.2. It 
is easier to evaluate the difference Vg R \r; ///) — Vc+l(?") accurately, using eq. (81), than to directly 
evaluate Vc+h(r): 



T/(% \ -it r \ 2Cp f A - [1 + iqr + \{iqrf] 

V s ( r ^f) - V c+L{r) = Im / dq * + const. 

TT J c qr 



'C 3 

(148) 

Here, we choose the second scheme for a ( yj(q). The integral path C 3 is shown in Fig. 16, e.g. 
q = k + ik 2 (k — fif) 2 for < k < [if. We solve the RG equation for as(q) along this path similarly 
to above. We may ignore the r-independent constant on the right- hand- side of eq. (148). Then, 
subtracting eq. (148) from V s (r; ///) computed in the second scheme, we obtain Vc+h(r). 
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